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ABSTRACT 


In  situ  hydrographic  data  and  remotely  sensed  sea  surface  temperature  data  were 
analyzed  to  describe  the  kinematics  and  dynamics  of  a  cyclonic  eddy  observed  off  Point 
Arena,  California,  in  May  1993.  The  hydrographic  data  were  first  objectively  analyzed 
using  multiquadric-biharmonic  interpolation  (MQ-B)  to  obtain  a  three-dimensional  field 
of  density.  This  was  then  used  in  a  primitive  equation  (PE)  model  with  a  digital  filter 
initialization  method  to  arrive  at  the  dynamically  balanced  three-dimensional  velocity 
field,  including  the  vertical  motion.  Several  aspects  related  to  geostrophic  adjustment,  to 
data  analysis  and  to  the  lack  of  exact  synopticity  in  the  observations  were  examined;  these 
include  the  effects  of  internal  tides,  eddy  rotation  and  MQ-B  parameters  sensitivity 
analysis. 

The  results  reveal  a  cyclonic  circulation  with  meanders  in  the  otherwise  nearly  circular 
flow,  and  maximum  horizontal  speeds  of  about  40-50  cm/s  at  the  surface.  The  meanders 
have  associated  patterns  of  radially  aligned  patches  of  sinking  and  rising  motion,  extending 
coherently  to  about  500  m  depth,  with  maximum  vertical  velocities  of  20-22  m/day 
between  100  and  150  m  depth.  Because  of  the  strong  horizontal  currents  and  the  short 
horizontal  length  scales  involved,  the  particle  vertical  excursions  are  very  small. 

The  vertical  motion  in  the  domain  was  also  diagnosed  from  the  analyzed  density  field 
using  the  quasigeostrophic  (QG)  omega  equation.  The  resulting  vertical  velocities  are  in 
good  agreement  with  those  from  the  initialized  PE  model  except  that  the  QG  values  were 
about  30  percent  stronger  in  magnitude.  These  discrepancies  were  shown  to  be  caused  by 
the  neglect  of  ageostrophic  advection  in  the  QG  system. 


The  views  expressed  in  this  dissertation  are  those  of  the  author  and  do  not 
reflect  the  official  policy  or  position  of  the  Department  of  Defense  or  the 
U.  S.  Government. 
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I.  INTRODUCTION 


The  large  scale  and  the  mesoscale  circulation  in  eastern  boundary  currents  (EBC) 
are  becoming  well  established  areas  of  active  research.  This  is  particularly  true  in  the 
coastal  zone,  where  the  interaction  between  the  eastern  boundary  currents  and  the  boundary 
is  strongest.  Consequently,  a  series  of  research  experiments  designed  to  study  the  dynamics 
of  these  regions  are  being  carried  out.  An  area  that  has  seen  a  great  deal  of  research  is  off 
the  west  coast  of  the  United  States,  mainly  the  Oregonian  and  Northern  Californian  coasts 
(California  Cooperative  Oceanic  Fisheries  Investigation,  Coastal  Ocean  Dynamics 
Experiment,  Coastal  Transition  Zone  experiments,  among  others). 

An  extensive  work  describing  the  circulation  of  the  California  Current  System, 
assembling  results  from  various  early  experiments,  is  found  in  Hickey  (1979).  This  study 
documents  the  existence  of  a  broad,  slow  southward  flow  offshore,  with  strong  coastal 
up  welling  in  selected  areas  during  periods  of  strong  northerly  wind  events.  The  study  also 
describes  an  undercurrent  flowing  poleward  and  the  reversal  of  the  inshore  surface  current, 
that  flows  poleward  when  the  wind  relaxes. 

The  early  interest  in  eastern  boundary  currents  was  indeed  linked  to  the  study  of 
upwelling,  an  important  phenomenon  that  takes  place  over  large  areas  of  western 
continental  margins  and  which  contributes  to  the  enhancement  of  the  nutrient  content  of 
coastal  waters  (e.g.  Huyer,  1983  and  Huyer,  1984).  It  was  later  recognized,  mainly  with  the 
availability  of  satellite  imagery  (Rienecker  and  Mooers,  1989;  Pares-Sierra  et  al. ,  1993), 
that  eastern  boundary  currents,  once  considered  slow  and  steady  flows,  are  in  fact  very  rich 
in  mesoscale  features  mainly  resulting  from  the  interaction  of  the  currents  with  the  coast, 
bottom,  and  local  winds.  From  this  interaction  there  results  an  intricate  circulation  where 
cold  filaments  extending  offshore  and  permanent  and  recurring  eddies  are  common 
features. 

With  the  identification  of  these  mesoscale  features,  the  ocean  sciences  community 
started  to  investigate  reasons  for  their  formation  and  development.  This  led  to  a  wealth  of 
observational  programs  utilizing  a  variety  of  remotely  sensed  and  in  situ  techniques,  not 
only  off  the  west  coast  of  the  United  States  but  also  in  other  eastern  boundary  current 
regions  (e.g.,  Ramp  et  al.,  1991;  Sousa,  1992;  Le  Groupe  Tourbillon,  1983).  The 
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importance  of  mesoscale  eddies  in  the  transport  of  heat,  salt  and  other  bio-physical 
properties  in  the  coastal  ocean  is  also  being  investigated. 

It  now  seems  clear  that  coastal  upwelling,  mesoscale  eddies  and  cold  filaments  are 
closely  related.  An  equatorward  coastal  jet  will  develop  meanders  by  virtue  of  instabilities 
in  the  flow,  associated  with  vorticity  and  curvature  mechanisms.  The  mean  relative 
vorticity  of  the  jet  increases  toward  the  coast,  and  the  ocean  depth  decreases  in  this 
direction,  assuming  a  sloping  bottom.  This  causes  the  mean  potential  vorticity  to  increase 
towards  the  coast,  so  the  meandering  of  the  jet  behaves  like  a  planetary  Rossby  wave, 
propagating  upstream  relative  to  the  flow.  This  propagation  speed  can  be  of  the  same  order 
of  magnitude  as  the  mean  equatorward  jet  speed,  and  thus  some  meanders  can  become 
stationary.  This  effect,  in  conjunction  with  local  topographic  features  such  as  capes,  which 
tend  to  induce  their  own  vorticity  effects  in  the  flow  and  also  cause  appreciable  local  wind 
changes,  is  a  good  scenario  for  the  appearance  of  permanent  or  recurrent  eddies  at  certain 
locations  along  the  coast.  These  mechanisms  may  produce  the  offshore  surface  flow  of 
recently  upwelled  waters,  that  tends  to  be  stronger  just  south  of  the  capes.  In  this  situation, 
an  eddy  could  entrain  and  advect  offshore  the  upwelled  coastal  waters. 

The  eddies  will  not  remain  stationary  once  generated.  According  to  the  analytical 
theory  of  Nof  (1981, 1984),  the  eddies  will  propagate  westward,  away  from  the  boundary. 
Although  Nof’s  analyses  are  based  on  somewhat  strong  assumptions,  the  general  results  are 
valid  in  a  qualitative  point  of  view  and  have  been  shown  to  be  in  good  agreement  with  some 
physical  cases  (Bowman,  1985;  Andrews  and  Scully-Power,  1976).  This  westward 
propagation  could  contribute  to  the  formation  of  the  filaments  observed  in  the  vicinity  of 
capes;  the  filaments  often  have  two  or  more  coupled  eddies  at  their  extremity. 

A  useful  technique  for  the  study  of  eddies  and  cold  filaments  is  numerical 
simulation  (e.g.  Haidvogel  et  ah,  1991;  Batteen  et  ah,  1992;  Pares-Sierra  et  al.,  1993).  In 
the  California  Current  region,  Paduan  and  Niiler  (1990)  and  Swenson  et  al.  (1992)  applied 
quasigeostrophic  theory  to  the  motion  of  clusters  of  surface  drifters.  Their  results  suggest 
that  meanders  in  the  seaward  flowing  jets  and  filaments  have  secondary  circulations  and 
vertical  motions  of  10-20  m/day  at  the  base  of  the  surface  mixed  layer.  Walstad  et  al. 
(1991)  used  a  quasigeostrophic  model  to  study  the  California  coastal  jet  observed  during 
the  1987  Coastal  Transition  Zone  (CTZ)  program.  The  focus  of  this  study  was  to  determine 
the  rotational  (nondivergent)  part  of  the  currents  including  the  “level  of  no  motion”  which 
was  found  to  be  shallower  toward  the  coast.  It  is  important  to  note  that  all  of  the  previous 
efforts  to  analyze  the  structure  and  dynamics  of  the  California  coastal  jet  and  eddies  have 
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been  based  on  quasigeostrophic  theory,  and  most  of  the  studies  such  as  those  based  on 
surface  drifters,  have  been  restricted  to  only  the  upper  mixed  layer. 

Other  authors  have  used  quasigeostrophic  theory  to  analyze  similar  situations  in 
other  ocean  areas.  In  the  study  of  vertical  motion  near  the  Gulf  Stream,  Lindstrom  and 
Watts  (1994)  used  the  quasigeostrophic  vorticity  equation,  among  other  methods,  to 
compute  vertical  motion  making  direct  use  of  observed  data,  since  their  data  had  sufficient 
temporal  resolution.  Henrick  et  al.  (1979)  made  use  of  an  analytical  model  based  on  the 
quasigeostrophic  potential  vorticity  balance  to  predict  perturbations  of  pressure,  density 
and  velocities  in  mesoscale  eddies.  Arhan  and  Colin  de  Verdiere  (1985)  used  data  collected 
during  the  Tourbillon  experiment  to  do  a  statistical  fit  of  streamfunction  and 
quasigeostrophic  vorticity  fields,  and  then  to  describe  the  dynamics  of  an  anticyclonic 
eddy.  The  vertical  motion  in  this  case  was  not  well  described.  A  number  of  other  studies 
present  only  a  description  of  data  collected  in  the  vicinity  of  eddies,  without  attempting  to 
go  beyond  basic  analyses. 

The  study  described  here  deals  with  such  an  eddy  observed  offshore  of  Pt.  Arena, 
California.  The  eddy  is  about  90  km  in  diameter  and  1000  m  depth.  The  eddy  was  observed 
between  22  and  25  May,  1993,  during  a  cruise  on  board  the  R/V  Point  Sur.  The  processes 
described  above  are  particularly  significant  in  the  area  where  the  eddy  was  observed;  the 
region  between  Pt.  Arena  and  Pt.  Reyes  has  the  most  favorable  upwelling  conditions  of  the 
entire  coast  (Huyer,  1984).  A  study  of  the  California  Current  System  off  Oregon  and 
Northern  California  by  Ikeda  and  Emery  (1984)  includes  satellite  imagery  showing 
filaments  and  eddies  offshore  of  Pt.  Arena  in  the  summer  of  1982,  suggesting  the  recurrent 
nature  of  the  eddy  observed  in  May  1993. 

The  objective  of  this  study  is  to  describe  the  complete  three-dimensional 
circulation  associated  with  a  closed  cyclonic  eddy  that  was  observed  in  the  California 
coastal  region.  The  description  is  accomplished  by  assimilating  a  unique  data  set  into  a 
primitive  equation  numerical  model.  The  methodology  used  is  to  objectively  analyze  the 
quasi-synoptic  hydrographic  data  collected  during  a  recent  Eastern  Boundary  Current 
(EBC)  research  program,  and  satellite  sea  surface  temperature  data,  and  then  to  use  the 
analyzed  density  field  as  input  to  the  model.  The  use  of  a  new  initialization  procedure  yields 
dynamically  balanced  velocity  and  mass  fields  that  are  used  in  the  description  of  the  eddy 
and  its  dynamics.  This  methodology  is  not  restrictive  to  the  application  above  (description 
of  the  observed  cyclonic  eddy)  or  to  the  California  coast.  It  can  be  applied  in  other  areas  of 
the  world  ocean  where  mesoscale  vertical  motion  has  been  poorly  or  not  determined  at  all, 
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as  long  as  the  available  data  has  enough  spatial  resolution.  These  areas  include,  for 
example,  other  eastern  boundary  current  regions  (e.g.,  the  Eastern  Atlantic  Current  along 
the  coast  of  Portugal).  Other  possible  applications  are  the  study  of  the  three-dimensional 
circulation  associated  with  coastal  jets  and  filaments,  and  frontal  zones. 

The  interpretation  of  the  results  from  the  model  includes  a  discussion  of  aliasing  or 
distortion  in  the  data  caused  by  processes  with  temporal  and  spatial  scales  of  the  same  order 
of  magnitude  of  the  spacing  of  the  hydrographic  stations.  Other  important  questions 
addressed  in  this  study  are  whether  internal  tides  were  present  during  the  collection  of  the 
hydrographic  data,  and  what  their  effect  might  be,  and  also  the  distortion  effects  related  to 
the  rotation  of  the  cyclonic  eddy  during  the  CTD  sampling. 

In  addition  to  presenting  a  complete  description  of  the  eddy,  the  hypothesis  that  the 
dynamics  of  the  eddy  can  be  interpreted  within  the  framework  of  quasigeostrophic  theory 
was  also  tested;  quasigeostrophic  theory  describes  adequately  the  dynamics  of  systems 
with  small  Rossby  numbers.  Attention  was  focused  on  the  vertical  motion,  for  several 
reasons.  First,  it  is  an  important  quantity  that  has  not  been  measured  directly  in  the  EBC 
program;  secondly,  it  is  responsible  for  the  vertical  advection  of  important  tracers  and  it 
influences  biological  productivity;  thirdly,  it  is  related  to  departures  from  geostrophy  and 
secondary  circulations;  finally,  it  is  essential  in  some  conceptual  models  of  vorticity  in 
meandering  jets  and  eddies. 

This  study  is  organized  as  follows:  Chapter  II  has  two  parts;  in  the  first  part  the 
available  data  is  described,  including  hydrographic  and  moored  data;  the  second  part 
presents  the  data  analysis,  including  the  interpolation  to  a  regular  grid  for  model 
applications.  Chapter  HI  describes  the  methods  used  for  the  computation  of  the  vertical 
velocity,  namely  the  circulation  model  initialized  in  such  a  way  as  to  isolate  the  slow  mode 
captured  by  the  data,  and  the  quasigeostrophic  omega  equation.  Chapter  IV  has  two  parts: 
the  first  part  presents  and  discusses  the  results  obtained  with  these  methods  and  describes 
the  physical  features  of  the  eddy.  It  also  includes  discussions  on  the  initialization  and 
geostrophic  adjustment,  and  potential  data  related  problems  such  as  the  effect  of  internal 
tides,  eddy  rotation  and  sensitivity  to  data  interpolation.  The  second  part  presents  and 
discusses  the  results  obtained  with  the  quasigeostrophic  omega  equation.  Finally, 
Chapter  V  gives  conclusions  and  final  recommendations. 
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II.  DATA 


A.  DATA  DESCRIPTION 

The  data  used  in  this  study  consists  of  satellite  data  and  hydrographic  data  from  a 
quasi-synoptic  Eastern  Boundary  Current  (EBC)  CTD  survey  off  Pt.  Arena,  California,  in 
May  1993.  The  objective  of  the  survey  was  to  observe  and  resolve  eddies  and  filaments  that 
had  been  seen  off  Pt.  Arena  in  earlier  studies  such  as  the  Coastal  Transition  Zone  (CTZ) 
program  in  1987-1988.  As  shown  below,  the  survey  was  successful  in  mapping  a  very  large 
fraction  of  a  closed  cyclonic  eddy  that  had  just  cut  off  and  separated  from  the  coast. 

1.  Satellite  imagery 

The  satellite  data  used  in  this  work  consists  of  images  obtained  at  the  Naval 
Postgraduate  School.  The  raw  images  contain  all  five  channels  of  the  Advanced  Very  High 
Resolution  Radiometer  (AVHRR)  flying  on  the  National  Oceanic  and  Atmospheric 
Administration  (NOAA)  series  of  Earth  observing  satellites.  This  sensor  observes  the 
surface  of  the  Earth  in  the  visible  and  infrared  (IR)  regions,  using  wavelengths  between  3.7 
and  11.5  pm.  The  sea  surface  temperatures  (SST)  were  computed  by  means  of  a 
multispectral  algorithm  that  converts  brightness  temperatures  to  sea  surface  temperature 
and  corrects  for  the  effect  of  the  atmospheric  water  vapor  on  the  observed  radiances. 

The  images  were  processed  using  software  developed  at  the  University  of  Miami. 
This  processing  consisted  of  correcting  the  images  for  appropriate  navigation 
(geographical  correction)  and  of  remapping  the  entire  image  into  the  region  of  interest, 
centered  at  a  point  off  Point  Arena,  California.  The  ground  resolution  depends  on  the 
distance  between  each  particular  pixel  and  the  satellite  nadir  and  is  estimated  in  the  range 
of  1  to  3  km. 

Most  of  the  images  collected  during  the  period  of  the  cruise  contained  large 
amounts  of  cloud  cover;  for  this  reason,  only  the  images  for  21  and  22  May  1993  (23:00 
Greenwich  Mean  Time)  were  suitable  to  be  used  in  this  work.  An  example  of  an  image 
from  22  May,  obtained  from  only  the  channel  four  of  the  AVHRR,  containing  also 
information  on  the  area  covered  by  the  model  domain  (see  Chapter  HI)  and  hydrographic 
stations  (see  section  n.A.2),  is  given  in  Figure  1. 
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Figure  1.  National  Oceanic  and  Atmospheric  Administration  (NOAA)  satellite  image, 
obtained  from  Channel  4  of  the  Advanced  Very  High  Resolution  Radiometer 
(AVHRR),  of  the  Sea  Surface  Temperature  (SST)  off  Point  Arena,  on  22  May  1993, 
23  GMT.  The  crosses  represent  the  positions  of  the  CTD  stations;  also  shown  is  the 
area  of  study  (model  domain).  Temperatures  are  in  degrees  Celsius. 
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A  surface  temperature  map  was  obtained  by  sampling  the  multispectral  image  of 
May  22,  23:00  GMT  at  the  same  horizontal  positions  of  the  grid  used  in  the  model,  thus 
providing  a  grid  of  satellite-derived  SST  used  later  in  the  data  interpolation  (see 
section  n.B.a). 

2.  Hydrographic  data 

The  hydrographic  data  used  in  this  study  were  collected  during  a  Naval 
Postgraduate  School  cruise  on  board  R/V  Point  Sur,  between  21  and  25  May,  1993,  using 
a  Neil  Brown  Mark  El  Conductivity,  Temperature  and  Depth  (CTD)  instrument.  The  CTD 
casts  were  performed  at  an  array  of  stations  listed  in  Table  1.  The  column  labeled  MAX  Z 
means  the  maximum  depth  attained  in  the  cast;  the  column  labeled  DEPTH  is  the  depth  to 
the  bottom  at  the  station.  See  also  Figures  1  and  2. 


LONG  (W) 

LAT(N) 

StaNo 

DATE 

TIME  (h) 

MAXZ 

(m) 

DEPTH 

(m) 

-123.697 

38.546 

1 

22 

23.9 

479 

500 

-123.779 

38.643 

2 

22 

22.7 

419 

460 

-123.857 

38.750 

3 

22 

21.4 

371 

410 

-123.951 

38.859 

4 

22 

20.0 

415 

455 

-124.027 

38.970 

5 

22 

18.5 

471 

530 

-124.176 

38.889 

6 

23 

10.6 

2253 

2420 

-124.080 

38.782 

7 

23 

8.3 

2041 

2110 

-123.984 

38.675 

8 

23 

6.0 

1959 

1965 

-123.917 

38.587 

9 

23 

4.0 

1933 

1940 

-123.811 

38.478 

10 

23 

1.8 

1623 

1630 

-123.940 

38.393 

11 

23 

24.0 

2191 

2280 

-124.026 

38.493 

12 

23 

21.2 

3111 

3125 

-124.100 

38.591 

13 

23 

18.8 

2443 

2500 

-124.212 

38.700 

14 

23 

16.1 

3197 

3200 

-124.313 

38.808 

15 

23 

12.8 

3257 

3260 

Table  1.  Characteristics  of  the  CTD  Stations 
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LONG  (W) 

LAT(N) 

StaNo 

DATE 

TIME  (h) 

mmsm 

| 

-124.456 

38.703 

16 

24 

11.4 

2041 

3500 

-124.332 

38.595 

17 

24 

9.1 

2035 

3445 

-124.239 

38.507 

18 

24 

6.4 

3411 

ggjgjPI 

-124.131 

38.414 

19 

24 

mm 

2033 

3430 

-124.027 

38.318 

20 

24 

2.0 

2031 

3430 

-124.149 

38.236 

21 

25 

11.6 

2031 

3600 

-124.248 

38.338 

22 

25 

9.4 

2035 

3670 

-124.356 

38.422 

23 

25 

0.6 

3563 

3650 

-124.467 

38.525 

24 

24 

22.0 

2029 

3605 

-124.574 

38.621 

25 

24  . 

19.6 

2037 

3660 

-124.787 

38.807 

26 

24 

14.5 

2035 

3610 

-124.681 

38.714 

27 

24 

16.9 

2033 

3650 

-124.498 

38.356 

28 

25 

3.8 

2033 

3760 

-124.629 

38.282 

29 

25 

6.0 

2033 

3840 

-124.049 

38.150 

30 

25 

13.6 

2033 

3610 

-123.948 

38.061 

31 

25 

15.7 

2031 

3510 

-123.842 

37.966 

32 

25 

17.9 

2049 

3480 

Table  1.  Characteristics  of  the  CTD  Stations  (Continued) 


The  instrument’s  pressure,  temperature  and  conductivity  sensors  were  calibrated 
prior  to  and  after  the  cruise.  A  General  Oceanics  rosette  sampler  was  used  to  collect  water 
samples  in  situ  for  conductivity  calibration  after  the  cruise.  The  CTD  sampling  rate  was  32 
Hz.  The  raw  data  were  collected  and  processed  using  a  software  package  developed  at 
EG&G  Marine  Instruments.  The  processing  done  at  this  level  consisted  of  removing  bad 
data  points  (outliers)  using  a  first  difference  criteria  to  detect  suspicious  data  points  and 
linear  interpolation  to  perform  the  correction.  The  data  were  then  averaged  to  every  2  dbar 
and  recorded  in  files  for  later  data  analysis. 

An  underway  data  acquisition  loop  was  active  during  the  cruise,  automatically 
recording  averages  of  temperature  and  salinity  (among  other  data)  at  a  depth  of  2  m  every 
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Longitude  (deg) 


Figure  2.  Position  and  number  of  CTD  stations.  Current  meter  moorings  are 
indicated  by  the  symbol  ©.  Also  shown  are  SAIL  stations  (symbol  *),  the  1000  m 
isobath  and  the  area  of  study  (model  domain). 


30  seconds  (Ship’s  Data  Acquisition  System,  SDAS).  These  surface  data  were  used  to  help 
fill  data  void  areas  in  the  interpolation  phase,  and  will  be  referred  to  hereafter  as  SAIL  data 
(from  Serial  Ascii  Interface  Loop).  SAIL  data  were  selected  at  four  positions  (SAIL 
stations)  separated  by  approximately  the  same  spacing  of  the  CTD  stations,  which  is  about 
14  km  (see  Figure  2). 

The  SDAS  also  recorded  wind  speed  and  wind  direction.  Figure  3  shows  the  wind 
conditions  during  the  cruise.  Although  there  was  considerable  wind  variability  during  the 
survey  period,  there  is  no  indication  that  it  was  in  a  scale,  in  space  or  time,  comparable  to 
the  observed  ocean  eddy.  Thus,  although  they  may  be  important  in  the  early  eddy  formation 
stages,  the  local  effects  of  wind  and  surface  heat  flux  forcing  are  neglected  in  this  study. 
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Day  in  MAY  93 


Figure  3.  Top:  wind  vector  (direction  in  degrees  true);  Bottom:  wind  speed  in  m/s 
during  May,  1993  off  Pt.  Arena.  The  vertical  dotted  lines  mark  the  beginning  and 
end  of  the  cruise. 

A  hull-mounted  ADCP  looking  downward  recorded  velocity  data  continuously 
following  the  track  of  the  ship,  using  accurate  navigational  information  from  the  Global 
Positioning  System  (GPS)  to  compute  the  absolute  velocity.  The  ADCP  is  able  to  measure 
the  three  components  of  the  current  using  an  array  of  four  acoustic  beams;  the  Doppler  shift 
induced  in  the  returning  echoes  by  the  moving  column  of  water  is  processed  to  compute  the 
absolute  three  dimensional  current  vector.  An  appropriate  time  gating  of  the  echo  yields  a 
profile  of  the  velocity  field  in  the  vertical  at  several  depth  bins  (for  more  details  on  the 
ADCP  instrument,  see  Acoustic  Doppler  Current  Profilers,  1989). 

In  addition  to  the  ADCP  current  observations,  moored  current  meters  (conventional 
rotor  meters)  were  also  deployed  at  station  numbers  3,  7,  8,  9,  13,  17,  18,  19  and  23  (see 
Figure  2).  These  instruments  were  deployed  at  100, 150, 300  and  600  m.  One  mooring  also 
had  an  instrument  at  1800  m.  Only  current  meter  data  from  stations  7,8,9  and  13  were 
available  for  this  study.  The  remaining  instruments  were  recovered  at  a  later  time  and  their 
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data  are  still  being  processed  elsewhere  (Scripps  Institution  of  Oceanography  and  Oregon 
State  University).  Table  2  lists  the  vertical  distribution  of  instruments  by  mooring, 
identified  by  CTD  station  number. 


STATION 

NUMBER 

DEPTH  (m) 

100 

150 

300 

600 

1800 

3 

X 

X 

7 

X 

X 

X 

X 

8 

X 

X 

X 

X 

X 

9 

X 

X 

X 

13 

X 

X 

X 

X 

Table  2.  Spatial  Distribution  of  Current  Meters 


It  is  assumed  that  the  uncertainty  associated  with  ADCP  velocity  measurements, 
which  can  be  as  large  as  10  to  15%  of  the  mean  horizontal  current  (e.g.,  Smith  and 
Morrison,  1989;  Geyer  and  Signell,  1990),  is  well  beyond  the  typical  magnitude  of  low 
frequency  vertical  motions,  thus  large  enough  to  disregard  the  use  of  ADCP  vertical 
velocity  estimates  in  this  study.  However,  the  observed  horizontal  velocities  from  both  the 
ADCP  and  the  current  meters  were  used  to  define  the  level  of  known  motion  (see  section 
m.A.2)  and  to  compare  model  results  with  observations  (see  section  IV.A.l).  Time  series 
of  temperature  and  velocity  data  collected  by  the  current  meters  were  also  used  in 
estimating  the  effect  of  internal  tides  on  the  results  (section  IV.A.3.a). 

B.  DATA  ANALYSIS 

An  inspection  of  the  maximum  depth  attained  on  each  CTD  station  in  Table  1 , 
“Characteristics  of  the  CTD  Stations”,  shows  that  whereas  the  deeper  stations  can  be 
referenced  to  a  common  depth  of  2000  m,  the  shallower  stations  (1  to  5,  8  to  10)  do  not 
reach  that  depth  due  to  the  shoaling  bottom.  For  reasons  that  will  become  clear  further  on 
(section  HI. A),  a  similarity  method  was  used  here  to  extend  those  stations’  profiles  down 
to  2000  m,  to  achieve  a  regular  three-dimensional,  flat  bottom  domain.  This  method 
searches  for  temperatures  in  the  nearest  complete  profile  equal  to  the  deepest  temperature 
in  the  profile  to  be  extended.  After  the  match  is  found,  the  lower  part  of  the  former  is  copied 
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to  the  latter.  The  similarity  method  is  attributed  to  Helland-Hansen,  but  no  citation  was 
found. 

The  now  complete  (common  depth)  CTD  data  were  then  analyzed  using  an  optimal 
interpolation  method  first  introduced  by  Hardy  (1971),  the  multiquadric-biharmonic 
(MQ-B)  interpolation  method,  that  has  seen  several  recent  applications  in  meteorology 
(Nuss  and  Titley,  1994)  and  oceanography  (Chumbinho,  1993).  Other  authors  have  applied 
this  method  to  other  fields  of  research,  including  the  one  for  which  the  method  was  initially 
designed,  the  representation  of  topography.  This  method  has  been  studied  by  several 
authors  (Hardy,  1990;  Sirayanone,  1988;  Nuss,  1994),  who  believe  that  it  presents 
advantages  over  other  common  interpolation  schemes  (e.g.,  Barnes  and  Cressman  objective 
analyses)  in  applications  where  irregularly  spaced,  scattered  data  must  be  mapped  into  a 
regular  grid.  The  method  behaves  very  benignly  in  data  void  areas,  that  is,  it  does  not 
introduce  unrealistic  scales  in  these  areas,  as  well  as  it  does  not  introduce  scales  smaller 
than  the  spacing  between  data  points.  Moreover,  the  simplest  form  of  MQ-B  interpolation 
(the  collocation  mode,  after  Hardy,  1990)  is  of  very  simple  implementation. 

The  purpose  of  using  this  method  in  the  data  analysis  was  to  combine  the  CTD  data 
below  the  surface  with  the  satellite  data  at  the  surface,  and  to  extend  the  observed  gradients 
horizontally  outside  the  CTD  area  to  fill  the  model  domain.  However,  preference  is  given 
to  the  horizontal  density  gradients  because  they  determine  the  vertical  shear  of  the  currents; 
this  limits  somewhat  the  propagation  of  surface  temperatures  to  deeper  levels,  which  is  also 
desirable.  The  satellite-derived  SST  are  representative  of  the  surface  skin  depth  only,  and 
thus  may  not  be  related  to  the  temperatures  at  deeper  levels.  Haney  et  al.  (1995)  show,  in  a 
study  that  attempted  to  recover  temperature  profiles  from  shallow  CTD  profiles  using 
empirical  orthogonal  functions,  that  CTD  casts  to  shallow  depths  less  than  100  m  are  not 
capable  of  representing  the  whole  water  column.  Thus,  the  main  result  of  this  MQ-B 
analysis  is  the  correct  representation  of  the  horizontal  density  gradients  determining  the 
vertical  shear. 

1.  Multiquadric-biharmonic  interpolation  method 

This  section  gives  a  description  of  the  MQ-B  method,  following  closely  that 
provided  by  Nuss  (1993).  Consider  a  three-dimensional  irregular  surface  containing  a  field 
sampled  at  N  observation  points.  For  each  observation  point  we  know  the  location,  (x,y,z), 
and  the  value  of  the  field,  H.  This  irregular  surface  can  represent  any  scalar  (and,  in  more 
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recent  applications,  also  vector)  field.  The  purpose  is  to  reconstruct  the  irregular  surface 
based  on  the  observed  data. 

The  interpolation  of  the  field  value  in  any  arbitrary  point  X  is  based  on  a  weighted 
sum  of  radial  basis  functions: 

N 

HiX)  =  £a,.0(X-X(.),  (2.1) 

i  =  l 

where  a,  is  the  weight  of  the  kernel  function  Q  whose  argument  is  the  radial  distance 
between  the  field  point  X  and  the  observation  point  X,.  In  the  multiquadric-biharmonic 
interpolation  method,  the  kernel  function  is  a  quadric  function: 

e(X-Xf)  =  [  (X-  X,-)2  +  c2]1/2.  (2.2) 

The  constant  c  is  a  small  arbitrary  number  to  keep  the  function  from  vanishing  at 
the  observation  points.  For  the  three-dimensional  interpolation  case,  the  kernel  function  is 
therefore: 

1/2 

e(X-X;)  =  [(x-Xi)2+  (y-yi)2+  (z-zf  +  c2]  .  (2.3) 

To  solve  for  the  weights,  we  apply  Equation  2.1  to  all  N  observation  points, 
resulting  in  a  set  of  N  equations  with  N  unknowns: 

N 

h(Xj )  =  Y,aiQi(xpyj’zP-  (2-4) 

*•= i 

In  matrix  notation,  the  solution  is  given  by: 

a  =  Q-1  •  Hf  (2.5) 

where  Q  is  the  square  matrix  of  equation  coefficients  and  Hf  is  a  vector  of  the  observed 
field  values.  This  matrix  can  be  inverted  using  any  of  the  mathematical  routines  available 
for  solutions  of  linear  systems  of  equations  (e.g.,  Unpack,  1978,  or  Matlab,  1994). 

Once  the  weights  are  known,  the  interpolated  solution  to  any  point  in  the  field  (in 
usual  applications,  regularly  spaced  grid  points)  is  then  given  in  matrix  notation  by: 

Hg  =  Qg  •  3,  (2-6) 

where: 
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(2.7) 


1/2 

Qg  =  Wxg-xi)2+  c yg-yi)2+  (zg-zf+c2]  , 

(xg,yg,zg)  being  the  coordinates  of  each  desired  field  point.  The  matrix  on  the  left  hand  side 
of  Equation  2.7  has  as  many  rows  as  there  are  field  (grid)  points,  but  the  number  of 
columns  is  always  equal  to  the  number  of  observation  points  N.  In  other  words,  the 
solution  depends  only  on  the  number  of  observation  points  available  and  is  independent  of 
the  size  of  the  grid  to  be  interpolated.  The  two-dimensional  problem  is  just  a  particular 
case  of  the  more  general  three-dimensional  case,  with  all  z  coordinates  set  to  zero. 

The  condition  that  the  sum  of  the  weights  be  equal  to  zero  (Hardy,  1990)  is 
implemented  by  adding  an  extra  row  and  column  to  Q  as  follows: 


[i-i]  0 


(2.8) 


where  the  vectors  have  length  N,  and  by  appending  a  zero  to  the  data  vector  Hf.  This,  in 
turn,  introduces  a  new  term  a0  in  the  weight  vector: 


a2 


aN 


(2.9) 


The  new  term  a0  represents  the  average  value  of  the  entire  observation  data  and  thus 
must  be  added  to  the  final  interpolated  field.  The  multiquadric  matrix  in  Equation  2.6 
maintains  the  dimensions  mentioned  above  (i.e.,  MxN,  where  M  is  the  total  number  of  grid 
points),  since  for  the  field  interpolation  only  the  first  N  weights  are  used.  The  MQ-B 
interpolation  equations  are: 

Hg  =  Qg  •  a  +  ci0,  (2.10) 

with  i  -  1,...,  N  and  g  -  1,...,  M. 

The  method  presented  above  is  valid  for  noise-free  observations.  Unfortunately, 
real  data  include  observational  errors  that  must  be  taken  into  account.  This  is  achieved 
through  the  introduction  of  an  uncertainty  parameter  and  a  smoothing  parameter  (Nuss  and 
Titley,  1993)  in  Equation  2.4.  In  matrix  notation: 

Hf  =  [Q  +  N^ct2  •  8^.]  •  a  ,  (2.11) 
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where  X  is  the  smoothing  parameter,  a2  is  a  vector  of  uncertainties  in  the  observations 
(mean-squared  observation  error)  and  8y  is  the  Kronecker  delta  function.  The  difference 
now  is  that  the  diagonal  elements  of  Q  are  modified  to  include  a  factor  involving  those 
two  parameters.  Note  that  the  uncertainties  vector  can  have  different  values  for  different 
observation  points,  a  useful  property  that  gives  the  user  of  the  method  enough  freedom  to 
optimally  set  the  parameters  for  best  fitting  any  part  of  the  grid. 

2.  Data  interpolation 

The  multiquadric-biharmonic  interpolation  method  described  above  was  applied  to 
the  temperature  and  salinity  data  to  obtain  an  analyzed  three-dimensional  density  field  in 
the  model  domain.  Chumbinho  (1993)  presents  some  thoughts  about  the  MQ-B  method, 
namely  the  choice  of  the  interpolation  parameters  c  and  X,  and  the  choice  of  an  appropriate 
vertical  scale  exaggeration.  The  interpolation  parameters  were  set  to  values  close  to  those 
used  in  other  applications  of  MQ-B  (Nuss,  1994;  Nuss  and  Titley,  1994),  in  a  way  that 
guarantees  the  best  representation  of  local  maxima  and  minima,  without  distortion  of  the 
interpolated  field  in  regions  where  data  is  scarce  or  not  available.  The  vertical  exaggeration 
is  necessary  to  emphasize  the  correct  representation  of  the  horizontal  gradients  over  the 
vertical  ones,  which  are  usually  much  larger. 

The  number  of  temperature  observations  is  much  larger  than  the  number  of  salinity 
observations  because  of  the  use  of  satellite-derived  sea  surface  temperatures.  As  noted 
previously,  these  were  used  in  an  attempt  to  compensate  for  the  lack  of  CTD  data  in  several 
parts  of  the  domain.  For  this  reason,  the  data  interpolation  was  carried  out  differently  for 
these  two  state  variables. 

a.  Temperature  interpolation 

The  CTD  temperature  data  were  first  decimated  along  each  profile  to  the 
same  27  depth  levels  used  by  the  circulation  model  (-3.5,  -12.5,  -24.3,  -38.0,  -54.9,  -75.9, 
-101,  -131,  -166,  -206,  -252,  -303,  -360,  -423,  -493,  -570,  -653,  -744,  -842,  -948,  -1060,  - 
1182,  -1312,  -1450,  -1596,  -1750  and  -1916  m;  also,  see  section  IH.A).  This  decimation 
was  performed  following  the  filtering  criteria  defined  by  Crochiere  and  Rabiner  (1983)  for 
non-uniform  decimation,  and  according  to  the  idea  that  the  MQ-B  interpolation  works  best 
when  the  data  points  coincide  with  grid  points  (Hardy,  1990,  p.  164).  This  resulted  in  32 
in  situ  temperature  data  points  at  each  of  the  27  vertical  levels,  for  a  total  of  864  points. 
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The  satellite-derived  sea  surface  temperature  data  (SST)  were  sampled 
directly  from  the  multispectral  satellite  image  equivalent  to  the  single-channel  image  in 
Figure  1,  at  a  resolution  of  48x39  points,  respectively  in  the  alongshore  and  cross-shore 
directions,  covering  the  model  domain  represented  by  the  rectangle  in  Figure  2,  using  a 
geographical  to  Universal  Transverse  Mercator  (UTM)  coordinates  transformation  and 
posterior  rotation.  The  origin  of  the  UTM  coordinates  is  the  northwestern  comer  (or  lower 
left  comer,  in  the  figures).  The  total  number  of  SST  data  points  is  thus  1872,  which  is  much 
larger  than  the  32  temperature  data  points  at  the  second  level. 

The  MQ-B  performs  best  when  the  data  are  irregularly  spaced,  without 
sharp  contrasts  in  spatial  data  density  (Hardy,  1990).  This  was  readily  observed  in  the  first 
few  interpolation  runs  without  any  pre-processing:  the  three-dimensionally  interpolated 
temperature  field  presented  severe  distortions  in  the  upper  levels,  caused  by  the  contrast  in 
the  number  of  total  data  points  at  the  surface  and  at  the  level  immediately  below.  This 
problem  lead  to  the  necessity  of  pre-processing  the  temperature  observations. 

The  SST  data  were  thus  two-dimensionally  decimated  by  a  factor  of  seven, 
resulting  in  a  rectangular  array  of  8x7  data  points  covering  the  same  area,  again  following 
the  criteria  defined  by  Crochiere  and  Rabiner  (1983);  at  the  same  time,  a  two-dimensional 
MQ-B  interpolation  of  temperature  was  carried  out  at  each  level,  using  very  small 
parameters  (c=lxl0‘4,  X  infinitesimally  small),  to  artificially  increase  the  number  of 
observations  at  depth  from  25  to  64  (see  below).  The  interpolation  parameters  mentioned 
above  introduced  no  smoothing  in  the  interpolated  field,  and  caused  the  resulting  field  to 
follow  the  data  points  exactly  (the  interpolation  results  were  checked  using  three- 
dimensional  visualization  of  the  data  points  against  the  interpolated  field).  Only  the 
decimated  temperature  profiles  corresponding  to  stations  1  to  25  were  used  in  this 
two-dimensional  interpolation,  and  the  resulting  field  was  mapped  into  a  square  8x8  grid 
about  the  same  size  as  the  horizontal  area  occupied  by  these  stations,  as  illustrated  in 
Figure  4. 

The  last  step  in  the  temperature  interpolation  was  to  combine  the  8x7 
decimated  SST,  the  8x8x27  two-dimensionally  interpolated  temperatures,  and  the  7x27 
points  from  the  remaining  seven  CTD  stations  (numbers  26  to  32)  in  a  three-dimensional 
MQ-B  interpolation,  using  c-2xl0'4,  X=lxl0'3  and  an  aspect  ratio  of  about  1/60  (achieved 
by  scaling  the  horizontal  dimensions  to  unity  and  leaving  the  vertical  dimension  unaltered). 
The  resulting  three-dimensionally  analyzed  temperature  field  followed  the  data  points 
closely;  profiles  of  analyzed  data  were  compared  with  nearby  CTD  profiles,  showing  a 
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Distance  Alongshore  (km) 


Figure  4.  Schematic  representation  of  the  2-D  temperature  interpolation  at  each 
level.  The  dots  are  the  model  domain  grid  points,  the  crosses  are  the  CTD  stations 
and  the  rectangle  is  the  area  covered  in  the  2-D  interpolation,  mapped  by  a  8x8 
grid. 

good  match.  The  maximum  differences  between  analyzed  and  observed  global  maximum 
and  minimum  temperatures  were  0.11  and  0.07  degrees  Celsius,  respectively.  This 
analyzed  field  had  48x39  points  at  each  of  the  27  levels. 

b.  Salinity  interpolation 

The  MQ-B  interpolation  of  salinity  was  considerably  simpler  than  the 
interpolation  of  temperature,  because  there  were  no  sharp  contrasts  in  the  number  of 
observations  between  levels.  The  CTD  salinity  profiles  were  non-uniformly  decimated  in 
the  vertical,  in  a  manner  similar  to  the  temperature  vertical  decimation,  to  the  same  27 
depth  levels,  resulting  in  32x27  observation  points.  Four  other  salinity  observation  points 
were  included  in  the  salinity  data  vector,  corresponding  to  the  SAIL  stations,  and  surface 


17 


salinities  were  taken  from  an  average  of  the  first  two  CTD  records  on  each  profile.  The  total 
number  of  salinity  observations  was,  therefore,  900  data  points. 

These  observations  were  three-dimensionally  MQ-B  interpolated  using 
c=lxl0'4,  ?i=lxlO-3  and  an  aspect  ratio  of  about  1/60,  similar  to  the  temperature 
interpolation.  The  slightly  different  value  for  c  yielded  better  results,  maybe  because  of  the 
different  spatial  distribution  of  salinity  observation  points  compared  to  that  of  the 
temperature  observations.  The  resulting  three-dimensionally  analyzed  salinity  field  also 
followed  the  data  points  closely,  with  very  good  match  between  analyzed  profiles  and 
nearby  CTD  salinity  profiles.  The  maximum  differences  between  analyzed  and  observed 
global  maximum  and  minimum  salinities  were  0.01  and  0.04  psu,  respectively.  The 
dimensions  of  this  field  were  the  same  of  the  temperature  field  (48x39x27). 

3.  Analyzed  three-dimensional  density  field 

Once  the  temperature  and  salinity  data,  from  the  various  sources  described  above, 
were  interpolated  and  analyzed,  the  three-dimensional  (48x39x27)  density  field  was 
computed  using  the  formulae  for  p(S,T,p)  in  the  International  Equation  of  State  (IES  80) 
found  in  Gill  (1982).  The  values  of  density  were  then  converted  to  sigma  values: 

a  (S,  T, p)  =  [p(S,  T, p)  -  1000]  (2.12) 

where  S  is  salinity,  T  is  temperature  and  p  is  pressure. 

The  result  was  a  three-dimensional,  analyzed  density  field  from  which  the  three- 
dimensional  currents  were  diagnosed  using  the  primitive  equation  model,  and  the 
quasigeostrophic  vertical  motion  was  computed  from  the  quasigeostrophic  omega 
equation.  An  implicit  assumption  of  synopticity  in  the  data  was  made  throughout  the 
analysis.  The  hydrographic  data  were  collected  during  a  period  of  three  days,  during  which 
several  processes  can  act  to  distort  the  sampled  field.  Among  these  are  the  rotation  and 
translation  of  the  nearly  circular  eddy  (Figure  1),  and  the  propagation  of  internal  waves  in 
the  domain  sampled  by  the  CTD.  The  most  important  internal  waves  that  could  be  present 
are  internal  tides,  which  have  been  known  to  contaminate  data  sets  in  other  areas.  The 
possible  effects  of  these  processes  are  analyzed  and  discussed  in  section  IV.A.3.a. 
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III.  METHODS 


After  the  data  were  analyzed  and  the  three-dimensional  density  field  computed,  a 
field  of  dynamically  balanced  currents  in  the  area  of  study  was  obtained  by  initializing  a 
primitive  equation  model.  To  help  interpret  the  results,  the  vertical  velocity  was  also 
diagnosed  from  the  (quasigeostrophic)  omega  equation.  Both  the  primitive  equation  model 
and  the  quasigeostrophic  omega  equation  were  applied  in  the  volume  defined  by  the 
analyzed  density  field  (48x39x27,  respectively  alongshore,  cross-shore  and  in  the  vertical) 
whose  horizontal  extension  is  depicted  in  Figure  1  (hereafter  domain).  This  domain  was 
chosen  to  cover  a  larger  horizontal  area  than  the  domain  defined  by  the  CTD  stations  to 
allow  for  model  boundary  effects  and  possible  eddy  translation. 

This  chapter  describes  several  aspects  relevant  to  each  of  these  solution  methods. 

A.  CIRCULATION  MODEL 

The  objectively  analyzed  density  field  was  used  in  a  circulation  model  for  the 
diagnostic  computation  and  analysis  of  the  dynamically  balanced,  three-dimensional 
velocity  and  density  fields  corresponding  to  the  slow  mode  assumed  captured  in  the 
observations. 

1.  Description 

The  model  used  in  this  study  is  the  multi-level,  primitive  equation  model  presented 
by  Haney  (1974);  it  utilizes  a  B  grid  in  the  horizontal  and  sigma  coordinates  in  the  vertical, 
with  vertical  staggering  of  the  vertical  component  of  velocity.  The  inviscid,  hydrostatic  and 
Boussinesq  assumptions  are  used  to  simplify  the  set  of  governing  equations.  This  means 
that  the  ocean  is  considered  incompressible,  and  density  is  replaced  by  a  constant  except 
where  it  is  multiplied  by  the  acceleration  of  gravity. 

This  model  has  been  improved  with  the  modifications  introduced  in  the  convective 
adjustment  algorithm  by  Adamec  et  al.  (1981)  and  in  the  numerical  computation  of  the 
pressure  gradient  force  in  sigma  coordinates  and  in  the  vertical  discretization  scheme  by 
Haney  (1991).  The  convective  adjustment  scheme  imposes  dynamic  stability  by  using  an 
algorithm  that  conserves  heat  and  momentum  and  mixes  these  two  properties  between 
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levels  at  the  same  ratio,  requiring  a  final  Richardson  number  (after  the  scheme  has  been 
applied)  greater  than  the  critical  value  of  0.25.  The  new  vertical  discretization  and  pressure 
gradient  force  computation  were  designed  to  improve  this  computation  over  steep 
topography. 

In  this  study,  the  dynamics  of  the  observed  eddy  will  be  studied  under  the  working 
hypothesis  that  surface  forcing  (wind  stress  and  buoyancy  fluxes)  and  variations  in  bottom 
topography  can  be  neglected.  This  simplifies  the  model  by  permitting  the  use  of  a  single 
thermodynamic  variable,  o,  given  in  Equation  2.12,  and  a  flat  bottom.  The  possible  effects 
of  these  assumptions  will  be  discussed  in  Chapter  V. 

2.  Initial  conditions 

To  achieve  dynamic  initialization  of  the  primitive  equation  model,  information  on 
the  initial  mass  and  velocity  fields  is  needed.  The  former  is  the  objectively  analyzed 
density,  or  o;  the  latter  can  be  set  in  one  of  two  ways,  if  no  direct  current  observations  are 
available  in  the  model  domain.  The  initial  shear  velocity  field  can  either  be  set  to  zero, 
implying  that  no  a  priori  knowledge  of  the  initial  velocities  is  available,  or  it  can  be 
specified  from  the  initial  mass  field  using  the  hydrostatic  and  geostrophic  relations.  These 
two  initial  velocity  fields  have  different  implications  in  the  final  balanced  slow  mode 
solution,  as  will  be  seen  in  section  IV.A.2. 

The  depth  averaged  flow  is  also  computed  from  the  hydrostatic  and  geostrophic 
relations,  assuming  a  level  of  known  motion  determined  by  the  current  meter  data.  The 
resulting  vertically  averaged  flow,  being  nondivergent,  is  then  used  to  compute  a 
streamfunction  for  the  vertically  averaged  velocities.  This  streamfunction  is  then  kept 
constant  during  the  short  model  initialization  run. 

The  level  of  known  motion  was  chosen  to  be  2000  m,  with  a  current  of  2  cm/s 
directed  alongshore  to  the  northwest,  based  on  current  meter  data.  This  apparently  simple 
method  was  preferred  over  a  more  complex  method  of  determining  a  level  of  known  motion 
by  using  the  horizontal  ADCP  currents.  This  approach  would  have  required  an  objective 
analysis  of  the  ADCP  currents  and  the  solution  of  a  Poisson  equation  for  the  rotational  part 
of  the  velocity  field  with  consequent  specification  of  boundary  conditions  for  the 
streamfunction.  The  correct  determination  of  the  depth  averaged  flow  is  important  because 
of  its  effect  in  the  advection  of  relative  vorticity  perturbations;  however,  it  will  be  shown 
that  the  solution  is  barely  sensitive  to  variations  in  the  setting  of  the  bottom  current. 
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3.  Boundary  conditions 


The  model  has  open  lateral  boundaries  governed  by  modified  Ross  and  Orlanski 
(1982)  radiation  conditions,  hereafter  RO  boundary  conditions.  This  modified  scheme 
serves  to  radiate  away  from  the  domain  any  outgoing  waves  (mainly  gravity  waves),  and  to 
include  inward  advection  of  properties  specified  at  the  boundaries.  The  RO  scheme  can  be 
summarized  in  three  steps: 

•  computation  of  a  normalized  phase  velocity  8  for  propagating  waves  (see 
Equation  3.5  below); 

•  for  inflow  points  (8  <  0),  the  vorticity,  divergence  and  density  are  set  to  the 
observed  boundary  values; 

•  for  outflow  points  (8  >  0),  the  vorticity,  divergence  and  density  are  advected  from 
the  interior  using  equations  of  the  type: 


Bt  Bn 


(3.1) 


In  Equation  3.1,  \  is  either  divergence,  vorticity  or  density,  and  the  spatial 
derivative  is  taken  along  an  outward  normal  n  to  the  boundary.  Once  the  values  of 
divergence  and  vorticity  are  known  at  the  boundaries,  their  definitions  are  used  to  compute 
the  normal  and  tangential  components  of  the  baroclinic  horizontal  velocity.  For  outflow, 
density  is  also  advected  along  the  boundary  with  an  equation  similar  to  Equation  3.1,  with 
the  spatial  derivative  generalized  to  include  its  tangential  counterpart,  and  a  new  phase 
velocity  82  computed  as  the  spatial  four-way  average  of  the  surrounding  tangential 
velocities  at  each  boundary  point. 

The  boundary  equations  for  vorticity  and  divergence  in  the  outflow  case  were 
modified  from  the  original  equations  given  in  RO,  since  their  equations  were  mainly 
designed  to  deal  with  gravity  waves,  with  divergence  waves  propagating  faster  than 
vorticity  waves.  RO  use  the  phase  speed  8  to  advect  vorticity,  and  Ax/ At  to  advect 
divergence,  where  Ax  is  the  grid  spacing  and  At  is  the  time  increment.  This  means  that  the 
divergence  at  the  boundary  is  set  equal  to  the  value  at  the  adjacent  interior  point,  computed 
from  the  prognostic  equations  governing  the  interior. 

The  velocity  field  in  the  model  initialization  runs  described  below  presents 
noticeable  inertial  oscillations,  where  divergence  and  vorticity  waves  propagate  at 
comparable  speeds.  The  different  advection  speeds  for  one  and  the  other  in  the  RO  method 
caused  the  components  of  horizontal  velocity  to  become  excessively  large  on  the 
boundaries.  Because  of  this  it  was  necessary  to  include  the  inertial  terms  and  a  damping 
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The  computation  of  an  appropriate  phase  speed  8  was  also  revised  from  RO. 
Following  the  work  of  Durran  et  al.  (1993),  the  computation  of  the  phase  speed  was 
changed  to: 

8  =  Ca  +  V  •  n  (3.5) 


where  V  is  the  velocity  vector  at  the  boundary,  n  is  an  outward  unit  vector,  and  Ca  is  a 
constant,  set  as  to  guarantee  propagation  of  waves  out  of  the  domain  during  most  of  the 
integration  time  (equivalent  to  the  shallow  water  phase  speed  of  gravity  waves).  V  •  n  is  a 
four-way  average  of  the  normal  velocity;  for  instance,  for  the  western  boundary  V  •  n  is 
computed  as: 


V  *  n  =  +  4j- 1  +  ulj  +  « u  - 1 1 


The  superscript  (n)  means  values  of  the  variables  at  time  n.  This  makes  the 
advection-type  Equation  3.1  a  completely  forward-in-time  scheme  (as  it  is),  which  is 
acceptable  for  upstream  advections  (it  is  used  to  compute  unknown  points  at  the  boundary 
from  known  interior  points,  in  outflow  situations).  The  final  vorticity  and  divergence 
equations  are,  respectively  (in  finite  differences  form): 


c+1>  -  C’  td"?,  -ti"’  ]  -/A<4"+1|-V|A!1C<") 
£>1"+1)  =Din)  -8^7  [fli"?i  -£>!”’]  +/A,^"+,)-7  |A^ln) 
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where  the  superscripts  (n)  and  (n+1 )  mean  time  level  n  and  time  level  n+1,  respectively, 
and  the  subscripts  B  and  B+l  mean  the  boundary  point  and  the  adjacent  interior  point. 

Note  that  because  Equation  3.7  uses  implicit  time  differencing  for  the  Coriolis 
terms,  it  has  to  be  solved  simultaneously  for  the  new  time  levels.  This  is  done  at  each  time 
step  along  all  four  boundaries,  except  in  the  comers.  The  vorticity  and  divergence  in  the 
comers  are  set  equal  to  the  average  of  the  two  adjacent  points  that  have  already  been 
computed  by  advection  normal  to  their  boundaries. 

For  use  in  the  case  of  inflow  at  the  lateral  boundaries,  the  vorticity  was  specified 
from  its  initial  value  and  the  divergence  (also  from  its  initial  value)  was  zero.  These 
conditions  were  held  constant  during  the  short  model  initialization  runs. 

The  boundary  conditions  for  the  vertical  velocity  in  this  study  were  zero  at  the 
bottom  and  at  the  surface. 

4.  Digital  filter  initialization 

The  objectively  analyzed  density  field  and  geostrophically  balanced  currents  were 
used  as  initial  conditions  in  the  model.  Separate  initialization  runs  were  also  made  with  the 
initial  currents  at  rest.  These  fields  are  not  in  balance,  or  at  least  not  in  the  state  of  near 
balance  that  characterizes  a  quasigeostrophic  ocean.  A  digital  filter  initialization  (DFI) 
procedure  introduced  by  Lynch  and  Huang  (1992)  was  used  to  arrive  at  a  dynamically 
balanced  field  of  density  and  currents.  The  DFI  procedure,  which  removes  high-frequency 
oscillations  from  the  model  fields,  consists  of  performing  two  short  (12  h),  adiabatic 
(physics  turned  off),  forward  and  backward  integrations  of  the  equations,  starting  with  the 
same  initial  conditions.  The  resulting  model  variables  at  each  time  step  from  t  =-12  h  to 
t  =12  h  are  multiplied  by  a  spatially  uniform,  time  varying  weight.  The  time  series  of 
weights  corresponds  to  the  coefficients  of  a  digital  lowpass  filter  effectively  being  applied 
to  the  computed  fields,  i.e.: 

N 

q(n)  =  £ q(k)h(n-k )  (3.8) 

k  =  -N 

applied  only  for  n  equal  to  zero.  In  the  equation  above,  q  is  any  of  the  variables  w,v  or  o, 
while  h(n)  is  the  sequence  of  filter  coefficients,  and  N  is  the  number  of  time  steps  (equal  to 
the  integration  time  forward  or  backward  divided  by  the  time  step,  which  was  300 
seconds).  This  procedure  is  illustrated  in  Figure  5.  The  figure  shows  an  example  of  a  time 
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Figure  5.  Example  of  the  DFI  procedure  applied  to  a  typical  vertical  velocity  time 
series.  The  solid  line  is  vertical  velocity  in  m/day;  the  dash-dot  line  is  the  sequence 
of  filter  coefficients  (scale  on  the  right)  multiplying  that  time  series. 


series  of  q ,  in  this  case  the  vertical  velocity,  at  one  of  the  grid  points,  and  the  time  filter 
that  is  to  be  applied  to  it.  The  output  from  Equation  3.8  is  q  (0) ,  the  filtered  vertical 
velocity  at  time  zero  at  that  grid  point. 

The  digital  filter  utilized  by  Lynch  and  Huang  (1992)  consists  of  a  simple  lowpass 
prototype  multiplied  by  a  Lanczos  window  to  reduce  truncation  effects.  The  same  authors 
have  already  proposed  other  filter  design  techniques  (Huang  and  Lynch,  1993).  The 
possibilities  in  this  area  are  broad,  since  various  filter  responses  and  designs  can  be  tested, 
as  demonstrated  by  Viudez  et  al.  (1994).  In  the  present  study,  the  filter  was  implemented 
as  in  Lynch  and  Huang  (1992): 


h{n)  = 


sin  [nn/  (N+  1)  ]  sin(n0c) 
nn/(N+l)  rnz 


(3.9) 
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where  0C  is  the  cutoff  digital  frequency.  When  this  filter  is  used  with  a  cutoff  period  equal 
to  the  total  time  spanned  by  the  two  integrations  (24  hours)  and  a  time  step  of  300  seconds, 
the  magnitude  of  the  frequency  response  shown  in  Figure  6  is  obtained.  It  is  seen  that  the 


Figure  6.  Frequency  response  of  a  lowpass  filter  with  cutoff  period  24  h,  sampling 
period  300  s  (Lanczos  window).  The  dotted  line  marks  the  digital  cutoff 
frequency.  The  dash-dot  line  marks  the  response  for  a  period  of  12  h. 


cutoff  period  is  damped  by  about  5  dB.  It  is  possible  to  design  more  efficient  filters  that 
would  achieve  better  frequency  responses  with  fewer  coefficients.  As  explained  above  in 
Equation  3.8,  the  number  of  coefficients  in  the  filter  is  directly  related  to  the  length  of  the 
forward  and  backward  integrations.  Thus,  fewer  coefficients  means  less  computer  time  in 
the  initialization  process.  On  the  other  hand,  it  also  means  less  time  for  ageostrophic  spin- 
up  (development  of  secondary  circulations  from  the  unbalanced  mass  and  velocity  fields; 
see  section  IV.A.2).  A  more  efficient  filter  design  (for  example,  using  Butterworth  or 
elliptic  lowpass  filters)  would  have  to  take  this  constraint  into  account  in  the  computation 
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of  the  coefficients  to  allow  enough  integration  time  for  spin-up  (usually  no  less  than  3  to  6 
hours,  setting  a  minimum  limit  of  36  to  72  coefficients  for  a  finite  impulse  response  filter). 

The  application  of  the  filter  to  the  model  output  does  not  need  to  be  done  within  the 
circulation  model,  as  long  as  enough  computer  storage  space  is  available  to  save  the  model 
output  at  each  time  step.  In  addition  to  using  different  filters,  the  idea  of  using  the 
dynamically  balanced  (filtered)  mass  and  velocity  fields  as  initial  conditions  for  a  second 
DFT  (integration-filter)  procedure  is  being  tested  and  verified  with  model  generated  data,  in 
a  separate  study. 

In  the  original  paper,  Lynch  and  Huang  (1992)  show  that  the  DFI  procedure  is 
equivalent  to  a  nonlinear  normal  mode  initialization,  with  the  advantage  of  having  a  much 
simpler  implementation.  The  normal  mode  initialization  compares  the  magnitude  of  the 
modes  resolvable  by  a  model  through  long-term  integrations  (gravity  wave  modes  versus 
Rossby  modes);  the  same  modes  are  fitted  to  the  analyzed  data,  and  the  components  that 
have  negligible  contributions  to  the  long-term  solution  (i.e.,  the  gravity  waves)  are 
removed  from  the  data.  The  nonlinear  variant  of  the  normal  mode  initialization  takes  care 
of  the  re-generation  of  high  frequencies  that  takes  place  even  after  those  frequencies  have 
been  removed  from  the  data,  because  of  the  nonlinear  terms  in  the  model  equations. 
Therefore,  the  model  is  initialized  with  the  slowly  varying,  almost  nondivergent  basic  state 
captured  in  the  data.  The  DFI  procedure  removes  the  high  frequency  modes  with  the 
application  of  the  filter.  The  lowpass  filtered  fields  then  represent  the  dynamically  balanced 
slow  mode  at  the  time  of  the  observations. 

Horizontal  and  vertical  diffusion  was  turned  off  to  enable  the  short  backward  and 
forward  integrations  centered  at  time  zero.  When  compared  to  advective  processes  in 
regions  of  strong  mean  flows,  such  diffusion  is  negligible,  as  noted  by  Schmitz  and  Vastano 
(1975)  in  a  numerical  study  of  Gulf  Stream  rings  over  time  scales  of  two  months. 

B.  QUASIGEOSTROPHIC  OMEGA  EQUATION 

The  second  method  of  dynamical  analysis,  which  was  used  to  diagnose  only  the 
vertical  velocity  from  the  data,  was  to  solve  the  quasigeostrophic  (QG)  omega  equation. 
The  assumptions  allowing  the  simplification  of  the  complete  set  of  equations  governing 
geophysical  fluid  dynamics  into  QG  equations  on  a /-plane  are  the  following: 

•  Hydrostatic  balance; 

•  Geostrophic  velocity  and  vorticity  are  used  everywhere  except  in  the  divergence 
term  in  the  vorticity  equation; 
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•  Adiabatic  flow; 

•  Neglect  of  the  vertical  advection  of  vorticity  and  the  tilting  of  vorticity; 

•  Relative  vorticity  much  less  than  planetary  vorticity  in  the  divergence  term; 

•  Planetary  vorticity  is  constant; 

•  Static  stability  is  a  function  of  z  (depth)  only. 

An  important  concept  in  QG  theory  is  that  vertical  motion  is  entirely  due  to  the 
ageostrophic  part  of  the  currents,  because  the  geostrophic  currents  are  nondivergent.  Thus, 
using  the  assumptions  above,  information  on  the  vertical  motion  is  contained  in  the 
vorticity  equation  relating  local  changes  to  vorticity  advection  by  the  geostrophic  flow: 

55  s  =  -vvcs-/„v.v  .  (3-10) 


where  C,g  is  the  vertical  component  of  the  geostrophic  relative  vorticity  defined  as: 


dx 


(3.11) 


(the  last  equality  follows  from  the  curl  of  the  geostrophic  balance),  g  being  the 
acceleration  of  gravity  and / the  planetary  vorticity  (2£2sin<|)). 

Equation  3.10  has  been  used  to  estimate  vertical  motion  after  substituting  the 
continuity  equation  for  incompressible  fluids 

V«V  =  (3.12) 

dz 

for  the  horizontal  divergence.  The  resulting  equation  for  w  requires  knowledge  of  the  local 
variation  of  vorticity  with  time,  which  is  not  always  available.  A  more  useful  equation 
without  any  time  derivatives  (diagnostic  equation)  can  be  derived  by  combining  Equation 
3.10  with  another  equation  relating  vertical  velocity  to  hydrostatic  pressure  perturbations, 
the  buoyancy  equation  (Gill,  1982,  page  529): 

PoAr2w  +  [|.  +  vs.v]|=0  (3.13) 


where  N2  is  the  buoyancy  frequency  defined  as: 


t2  _  8 


fb  dz’ 


where  p  is  pressure,  p  is  density,  p0  is  a  reference  density,  w  is  the  vertical  velocity,  and  \g 
is  the  geostrophic  velocity  vector.  The  Laplacian  and  del  operators  are  two-dimensional 


(horizontal).  Equation  3.13  is  the  nonlinear  QG  version  of  the  hydrostatic  buoyancy 
equation. 

With  the  help  of  the  hydrostatic  relation 


(3.15) 


the  buoyancy  equation  can  be  written  as: 

gji  =  p0 N2w-gVg  •  Vp 


(3.16) 


From  the  thermal  wind  relations 


d^8  =  AA £ 

dz  p  (/dy 

= 

dz  p(/dx 


(3.17) 


an  expression  for  the  vertical  derivative  of  geostrophic  relative  vorticity  can  be  derived: 


3C 

-l8  - 

dz  ~ 


— l_V2n 

■  Pot  p 


(3.18) 


Taking  the  z-derivative  of  Equation  3.10,  the  horizontal  Laplacian  of  Equation 
3.16,  using  the  results  in  Equations  3.12  and  3.18,  rearranging  terms  to  eliminate  the  time 
derivative  and  using  the  QG  assumptions  of  constant  planetary  vorticity  (f0 )  and 
horizontally  constant  stability,  the  quasigeostrophic  equation  for  vertical  velocity  is 
obtained  (omega-equation): 

=  /,5(V'’V?,)  +  ^V!<v,,Vp)  ■  <3I9) 


or,  in  terms  of  the  geostrophic  vorticity  defined  in  Equation  3.11: 


Equation  3.20  shows  that  vertical  motion  can  be  diagnosed  (no  time  derivatives) 
from  the  instantaneous  pressure  field  only,  since  this  variable  determines  both  the 
geostrophic  currents  and  the  geostrophic  relative  vorticity  in  the  QG  framework.  The  term 
on  the  left  hand  side  contains  only  second  derivatives  of  vertical  velocity;  for  disturbances 
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with  wave-like  character,  it  is  proportional  to  -w.  The  contribution  of  each  of  the  terms  on 
the  right  hand  side  of  this  equation  to  the  vertical  motion  will  be  discussed  in  section  IV.B. 

The  data  available  for  this  study  is  a  snapshot  of  the  density  field,  with  considerable 
detail  in  both  the  vertical  and  the  horizontal  (after  the  analysis  described  in  section  II.B). 
This  is  an  ideal  situation  for  the  application  of  Equation  3.19  for  the  computation  of  vertical 
motion.  However,  it  presents  some  problems  because  the  forcing  terms  on  the  right  hand 
side  (differential  vorticity  advection  and  Laplacian  of  thickness  advection)  sometimes  tend 
to  cancel  each  other,  leading  to  errors  in  the  numerical  procedure  to  compute  the  response 
on  the  left  hand  side.  A  more  stable  alternative,  the  Q-vector  form,  was  introduced  by 
Hoskins  et  al.  (1978): 

(JV2V2+^JL)w  =  2V«Q,  (3.21) 

3  z 

where 


Q 


8_ 

Po 


r3Vs  „  ' 

a;  #vp 

L3y  * Vp 


(3.22) 


This  form  of  the  to-equation  (Equation  3.21  above  is  the  so-called  Trenberth  form, 
where  deformation  terms  arising  from  taking  the  Laplacian  of  a  product  have  been 
neglected)  shows  that  a  convergent  Q-vector  field  implies  upward  vertical  motion. 
Moreover,  the  vectors  themselves  are  parallel  to  the  ageostrophic  velocities  responsible  for 
the  vertical  motion.  This  Q-vector  form  presents  additional  advantages  since  the  forcing 
term  is  exact  within  the  QG  context.  Neglecting  the  deformation  terms  could  be  a  source  of 
error  in  the  vicinity  of  strong  jets  or  frontal  systems  (Carlson,  1991),  however  that  is  not  a 
problem  in  the  present  case  because  even  a  coastal  jet  has  sufficiently  small  velocity 
magnitudes  in  an  eastern  boundary  current.  The  forcing  functions  can  also  be  evaluated  on 
a  single  level  surface. 

Several  aspects  of  the  QG  assumptions  can  cause  the  results  from  the  w-equation 
to  be  either  under-  or  over-estimated.  The  departures  from  QG  arise  in  regions  where: 
1)  diabatic  effects  are  strong;  2)  static  stability  is  not  constant  in  the  horizontal; 
and  3)  ageostrophic  motions  due  to  curvature  are  so  strong  that  neglecting  ageostrophic 
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advection  becomes  an  important  source  of  error.  The  relative  importance  of  these  terms  will 
be  considered  in  the  discussion  of  the  QG  results  (section  IV.B). 

The  o)-equation  was  solved  using  both  the  original  MQ-B  analyzed  density  field, 
and  the  dynamically  balanced  density  field  resulting  from  the  DFT  procedure  of  section 
m.A.4,  to  compute  the  forcing  (right  hand  side  of  Equation  3.21).  In  both  solutions  the 
condition  of  w=0  on  all  boundaries  was  used. 
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Figure  7.  Horizontal  cross  section  of  the  dynamically  balanced  horizontal  velocity 
field  at  a  depth  of  100  m,  in  cm/s  (red  arrows).  The  longest  arrow  is  about  50  cm/s. 
The  symbol  ®  represents  the  CTD  stations.  The  yellow  arrows  are  the  current  vectors 
from  the  current  meters  at  the  same  depth  (stations  7.  8,  9  and  13).  The  white  line  is  a 
streamline  of  the  nondivergent  component  of  the  velocity.  The  dashed  lines  are 
included  for  reference  and  represent  the  vertical  cross-sections  shown  in  other  figures. 
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The  mean  flow  is  not  circular  around  an  hypothetical  eddy  center;  instead,  it  has  meanders, 
with  at  least  five  well  defined  troughs.  These  smaller  scale  features  (hereafter  sub-eddy 
scales)  have  been  found  in  other  systems,  particularly  in  frontal  systems  (e.g.,  Fiekas  et  al. , 
1994)  where  submesoscale  circulations  are  part  of  the  mesoscale  meandering  of  the  front 
itself.  This  finding  is  partly  corroborated  by  an  inspection  of  Figure  1 .  This  figure  shows 
several  feeder  bands  of  water  being  entrained  into  the  eddy  circulation,  suggesting  that  the 
pattern  of  motion  around  the  eddy  is  not  completely  circular. 

A  comparison  was  made  between  the  DFI  solution  of  the  horizontal  velocities  and 
the  currents  observed  by  the  hull  mounted  ADCP  (section  II.A.2).  Figure  8  shows  a 
horizontal  map  of  the  dynamically  balanced  horizontal  currents  and  the  horizontal  ADCP 
vectors  for  the  bin  centered  at  100  m  depth.  In  general,  the  DFI  solution  and  the  ADCP 
currents  compare  well  in  magnitude  and  direction.  This  is  especially  so  in  the  northeast, 
northwest  and  southwest  parts  of  the  domain,  even  though  the  DFI  solution  in  the 
northeastern  part  of  the  domain  (upper  side  in  Figure  8)  presents  some  undesirable  effects 
from  the  flat  bottom  assumption  (see  below).  The  region  to  the  southeast  of  the  eddy  is  an 
area  where  practically  no  CTD  observations  were  made,  and  thus  the  solution  is  almost 
entirely  dominated  by  the  performance  of  the  MQ-B  method  in  data  void  areas,  which  was 
said  earlier  to  be  very  benign.  This  means  that  the  density  gradients  observed  by  the  CTD 
stations  at  the  southeast  edge  of  the  eddy  become  smaller  and  smaller  as  a  function  of 
distance  from  the  eddy,  and  the  general  circulation  will  therefore  tend  to  remain  weakly 
cyclonic.  An  exception  is  the  small  anticyclonic  feature  centered  at  (31,15),  which  is 
persistent  in  the  DFI  results  down  to  several  hundred  meters  and  was  never  sampled  by  the 
ADCP  (results  not  shown).  Since  this  position  coincides  with  a  CTD  station,  the  resulting 
anticyclonic  circulation  is  the  response  to  the  observed  density  gradients  between  stations 
22,  21  and  30,  extending  quite  deep. 

In  the  analysis  of  Figure  8  it  should  also  be  kept  in  mind  that  the  ADCP  observes 
the  total  velocity  with  a  measurement  uncertainty  that  is  highly  variable  and,  in  most  cases, 
of  the  same  order  of  magnitude  of  the  signal  that  is  being  computed  in  the  DFI  solution  (the 
ageostrophic  part  of  the  flow).  Also,  the  DFI  solution  is  representative  of  the  slow  mode  at 
the  time  of  the  observations,  whereas  the  ADCP  is  observing  all  frequencies.  The  ADCP 
data  plotted  in  Figure  8  was  spatially  averaged,  resulting  in  effective  time  averages  that 
range  from  15  minutes  to  about  2  hours. 

The  meanders  in  the  mean  flow  have  a  distinct  field  of  vertical  motion  associated 
with  them.  This  is  visible  in  Figure  9  which  shows  the  vertical  velocity  at  100  m  (colored 
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Figure  8.  Horizontal  cross  section  of  the  dynamically  balanced  horizontal  velocity 
field  at  a  depth  of  100  m,  in  cm/s  (red  arrows),  with  superimposed  ADCP  horizontal 
velocities  for  the  depth  bin  centered  at  100  m  (yellow  arrows). 
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Figure  9.  Horizontal  cross  section  of  the  dynamically  balanced  vertical  velocity  field 
at  a  depth  of  100  m,  in  m/day  (colored  field),  with  superimposed  horizontal  velocities 
at  the  same  level.  The  symbol  ©  represents  the  CTD  stations.  The  red  arrows  are  the 
current  vectors  from  the  current  meters  at  the  same  depth  (stations  7,  8,  9  and  13).  The 
black  bold  line  is  a  streamline  of  the  nondivergent  component  of  the  velocity.  The 
black  thin  lines  are  included  for  reference  and  represent  the  vertical  cross-sections 
shown  in  other  figures. 


field)  with  the  horizontal  velocity  (black  arrows)  and  the  same  streamline  as  in  Figure  7 
superimposed.  The  w  field  is  distributed  mostly  around  the  eddy  in  patches  of  radially 
aligned,  alternating  upwelling  and  downwelling  centers,  with  maximum  vertical  velocities 
at  this  level  of  about  ±  20  m/day.  The  average  vertical  velocity  within  a  radius  of  35  km 
from  the  center  of  the  eddy,  however,  is  in  the  order  of  0.1  m/day. 

These  sub-eddy  scale  patches  of  rising  and  sinking  motion  are  related  to  meanders 
in  the  otherwise  circular  eddy  and  can  be  explained  by  QG  theory  by  considering  the 
changes  of  vorticity  in  a  water  parcel  following  the  fluid  through  the  meanders  (Holton, 
1979).  Figure  10  shows  a  typical  meandering  system  with  cells  of  negative  and  positive 
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Figure  10.  Meandering  jet  with  alternating  patterns  of  negative  (regions  A  and  C) 
and  positive  (region  B)  vorticity.  The  changes  of  vorticity  going  from  A  to  B  to  C 
cause  downwelling  upstream  of  the  trough  and  upwelling  downstream  of  the  trough. 

vorticity.  In  region  A,  the  parcel  acquires  negative  relative  vorticity;  as  it  moves  to  region 
B,  it  encounters  regions  of  increasing  relative  vorticity.  The  QG  vorticity  equation 
(Equation  3.10)  can  be  written: 

f  =  -/V-V  =  /|.  (4.2) 

Hence,  from  A  to  B  the  left  hand  side  of  Equation  4.2  is  positive;  with  zero  boundary 
condition  at  the  surface,  the  right  hand  side  shows  that  w  is  negative  at  depth,  that  is,  from 
A  to  B  there  is  sinking  motion.  If  the  meander  is  not  growing  or  decaying,  the  zero  vertical 
velocity  point  will  be  exactly  at  the  trough,  because  in  this  case  the  left  hand  side  is  zero 
there. 


The  opposite  happens  as  the  parcel  moves  from  region  B  to  region  C;  the  change  in 
the  relative  vorticity  of  the  parcel  is  now  from  positive  values  to  negative  values,  causing 
the  left  hand  side  of  Equation  4.2  to  be  negative.  Using  the  same  boundary  condition  on  the 
right  hand  side  a  positive  value  for  the  vertical  velocity  is  now  found  at  depth.  Thus  an 
upwelling  cell  lies  between  region  B  and  region  C.  The  flow  is  thus  three-dimensional,  with 
the  particle  describing  not  only  the  meanders  in  the  mean  flow  but  also  some  vertical 
excursion  in  both  senses  (up  and  down).  In  this  simple  but  illustrative  model  the  vertical 
displacement  of  the  particle  is: 

®  inversely  proportional  to  the  magnitude  of  the  horizontal  mean  flow  between  A 
and  B  or  between  B  and  C; 

•  directly  proportional  to  the  magnitude  of  the  vertical  velocity. 

The  second  factor  is  proportional  to  the  rate  of  change  of  vorticity  following  the 
motion  and  the  depth  over  which  the  vorticity  change  occurs.  All  else  being  equal,  large 
vertical  scales  tend  to  have  stronger  vertical  velocities.  In  the  Gulf  Stream  meanders,  for 
instance,  vertical  motions  of  100  to  200  m/day  have  been  found  (Lindstrom  and  Watts, 
1994).  In  this  region  the  vertical  motion  may  be  enhanced  by  the  large  density  gradients 
across  the  meanders,  which  cause  the  changes  in  vorticity  to  take  place  over  appreciable 
changes  in  water  column  depth,  even  for  small  displacements  from  the  jet  axis.  The  rate  of 
change  of  vorticity  is  largest  in  regions  of  strong  changes  in  curvature,  i.e.,  for  smaller 
horizontal  scales. 

The  upwelling  and  downwelling  patterns  described  above  can  be  identified  in 
Figure  9.  The  patches  of  vertical  motion  are  correlated  with  the  meanders  of  the  mean  flow 
depicted  by  the  streamline,  with  downwelling  upstream  and  upwelling  downstream  of  the 
troughs.  Similar  patterns  of  vertical  motion  have  been  identified  in  the  meanders  of  the  Gulf 
Stream  (Leach,  1987;  Bower  and  Rossby,  1989;  Bower,  1989).  The  relative  magnitudes  of 
the  horizontal  and  vertical  velocities  also  show  that  the  vertical  excursion  of  the  water 
parcels  in  these  patches  is  small.  A  typical  horizontal  length  scale  for  the  patches  is  about 
20-25  km.  This  distance  is  covered  in  about  15  hours  at  a  speed  of  40  cm/s;  in  15  hours  a 
parcel  moving  at  this  speed  will  be  displaced  a  maximum  of  only  12  m  in  the  vertical. 

Another  important  aspect  of  the  description  of  this  eddy  is  its  vertical  extension. 
This  can  be  appreciated  in  Figures  1 1  and  12.  The  former  is  a  vertical  cross-section  of  the 
vertical  velocity  along  the  x-axis  (alongshore),  corresponding  to  the  horizontal  line  in 
Figure  7;  the  latter  is  a  cross-section  of  the  vertical  velocity  along  the  y-axis  (across  shore), 
corresponding  to  the  vertical  line  in  Figure  7.  Both  have  units  of  m/day,  although  the  color 
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Figure  11.  Vertical  cross  section  of  the  dynamically  balanced  vertical  velocity 
(colored  field)  in  the  alongshore  direction  (from  the  northwest  (left)  to  the  southeast 
(right)),  in  m/day,  with  depths  in  meters.  Also  shown  are  the  two-dimensional  (u,w) 
vectors  in  the  same  section  (black  arrows).  Note  that  the  vertical  axis  is  not  linear. 
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Figure  12.  Vertical  cross  section  of  the  dynamically  balanced  vertical  velocity 
(colored  field)  in  the  across  shore  direction  (offshore  to  the  left),  in  m/day,  with  depths 
in  meters.  Also  shown  are  the  two-dimensional  (v,w)  vectors  in  the  same  section  (black 
arrows).  Note  that  the  vertical  axis  is  not  linear. 


scale  is  not  the  same  in  the  two  figures;  also,  note  that  the  vertical  axis  is  not  linear  (the  top 
500  m  take  about  half  the  axis).  The  vertical  extent  of  the  eddy  is  shown  to  be  about  500  m 
(extending  in  some  points  to  almost  1000  m),  with  maximum  vertical  velocities  occurring 
between  50  and  150  m  depth.  The  patches  of  rising  and  sinking  motion  are  not  uniquely 
related  to  divergence  and  convergence  of  the  along-section  vectors  superimposed  on  the 
figures  because  there  may  be  compensation  in  the  across  section  flow.  However  it  can  be 
seen  that  the  strong  local  maxima  of  rising  motion  generally  have  associated  with  them 
divergence  of  the  velocity  vectors  above  and  at  the  same  level  while  the  local  minima  are 
generally  associated  with  convergence  above  and  at  the  same  level. 

It  is  also  noticeable  that  in  Figure  12  there  are  large  horizontal  vectors  in  the  far 
right  of  the  figure,  that  is,  the  boundary  of  the  model  closer  to  shore.  This  northeastern 
boundary,  and  the  northwestern  boundary,  are  regions  of  strong  tangential  velocities,  at 
least  in  the  upper  levels  as  can  be  seen  in  Figures  1  and  7.  This  causes  some  problems  at 
the  boundary,  which  use  an  advective  velocity  based  only  on  the  normal  flow  in  the 
boundary  condition  (section  III.A.3).  Furthermore,  this  is  also  the  region  where  the 
Helland-Hansen  method  was  used  to  extend  the  inshore  CTD  profiles  to  a  common 
maximum  depth  (section  II.B),  thus  introducing  horizontal  density  gradients  that  may  not 
be  real.  These  horizontal  gradients  cause  unrealistic  vertical  shear  in  the  horizontal 
velocities  through  the  thermal  wind  relations.  For  these  reasons,  the  model  results  in  the 
northeastern  part  of  the  domain  and  closer  to  the  shore  may  not  be  reliable,  and  the 
discussion  is  therefore  limited  to  regions  where  the  joint  effect  of  boundary  and  missing 
topography  is  clearly  nonexistent  or  minimized. 

2.  Geostrophic  adjustment 

It  was  mentioned  in  the  preceding  section  that  the  results  described  above  were 
obtained  with  the  DFI  procedure  in  a  model  whose  initial  baroclinic  velocity  field  is 
geostrophic.  Another  possibility  is  to  set  the  baroclinic  currents  to  zero  at  the  initial  time. 
This  different  initial  condition  will  also  yield  a  final  slow  mode;  the  choice  between  one 
and  the  other  is  related  to  the  process  of  geostrophic  adjustment  itself. 

Geostrophic  adjustment  is  the  process  by  which  an  initial  disturbance  in  a  rotating 
fluid  releases  part  of  its  potential  energy  into  kinetic  energy,  until  a  steady  equilibrium  is 
achieved  (Gill,  1982).  The  real  ocean  is  in  a  constant  process  of  geostrophic  adjustment  due 
to  multiple  forcing.  In  a  simulated  ocean  without  forcing,  the  geostrophic  adjustment  takes 
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place  during  the  spin-up  phase,  where  the  imbalance  between  the  initial  density  and  mass 
fields  is  geostrophically  adjusted.  During  the  process,  about  one  third  of  the  released 
potential  energy  is  converted  to  kinetic  energy  of  the  equilibrium  solution,  and  the 
remaining  two  thirds  are  dissipated  in  inertial  oscillations  around  the  equilibrium  state  or 
propagated  away  in  the  form  of  gravity  waves.  This  equilibrium  state  is  an  approximate 
geostrophic  balance,  that  is,  the  equilibrium  state  is  not  found  by  solving  the  steady-state 
geostrophic  equations.  For  a  rotating  eddy,  the  equilibrium  is  an  approximate  gradient 
balance  between  the  Coriolis  terms,  the  pressure  gradient  terms  and  the  acceleration 
(curvature)  terms. 

The  initial  disturbance  can  be  thought  of  as  containing  a  variety  of  length  scales;  the 
way  the  geostrophic  adjustment  affects  each  of  these  scales  is  determined  by  an  important 
horizontal  length  scale,  the  Rossby  radius  of  deformation  (Rossby  radius),  given  by: 
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where  c  is  the  phase  speed  of  the  disturbance  and  /is  the  Coriolis  parameter.  For  scales 
much  smaller  than  the  Rossby  radius,  rotation  effects  are  negligible  because  the  pressure 
gradient  is  much  larger  (short  waves  behave  as  gravity  waves  and  disperse  momentum 
rapidly  from  the  initial  disturbance),  whereas  for  scales  comparable  to  or  larger  than  the 
Rossby  radius,  rotation  effects  are  of  the  same  order  of  magnitude  of  the  pressure  gradient 
(long  waves  become  inertial  oscillations  and  propagate  much  slower).  The  time  scale  of 
the  adjustment  process  is  the  inertial  period  ( 2n/f=  19  h  at  38.6°)  at  a  distance  from  the 
initial  disturbance  comparable  to  the  Rossby  radius,  but  allowance  has  to  be  made  for  the 
oscillations  about  the  equilibrium  state. 

The  considerations  above  can  be  completed  with  an  analysis  of  the  non- 
dimensional  potential  vorticity  associated  with  the  initial  disturbance  (Gill,  1982): 
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where  T|  is  the  surface  elevation  and  H  is  the  water  column  depth. 

The  Rossby  radius  is  the  length  scale  at  which,  for  geostrophic  motion,  the  potential 
vorticity  terms  in  Equation  4.4  are  of  the  same  order  of  magnitude;  short  waves  are 
dominated  by  the  vorticity  term,  whereas  long  waves  are  dominated  by  the  pressure 
perturbation  term  (p ’  =  gpr|).  The  term  including  r|  is  related  to  changes  in  the  mass 
(pressure)  field;  the  term  including  C,  is  related  to  changes  in  the  velocity  field.  Therefore, 
for  an  initial  potential  vorticity  that  must  be  conserved  during  the  adjustment,  the  process 
takes  place  differently,  and  the  final  state  is  different,  depending  on  the  scale  of  the 
disturbance  compared  to  the  Rossby  radius.  For  scales  smaller  than  the  Rossby  radius, 
where  gravity  is  the  main  restoring  force,  the  final  relative  vorticity  dominates  the  potential 
vorticity,  and  the  mass  field  adjusts  to  be  in  geostrophic  balance  with  the  resulting  velocity 
field.  For  scales  comparable  to  or  larger  than  the  Rossby  radius,  where  rotation  is 
comparable  to  gravity,  the  final  mass  field  dominates  the  potential  vorticity,  and  it  is  the 
velocity  field  that  adjusts  itself  to  a  geostrophic  balance  with  the  mass  field.  This  is 

illustrated  in  Figure  13. 
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Figure  13.  Relation  between  the  mass  and  velocity  fields  during  the  geostrophic 
adjustment,  according  to  the  length  scale  (L)  relative  to  the  Rossby  radius  of 
deformation  (Rd).  V0  is  the  initial  velocity  field  and  p  is  the  initial  mass  field. 

Thus  it  is  seen  that  the  choice  of  the  initial  velocity  field  is  an  important  issue  for 
scales  smaller  than  the  Rossby  radius.  The  density  field  is  not  the  proper  information  with 
which  to  describe  these  scales,  because  it  adjusts  to  the  velocity  field.  Observations  of  the 
velocity  field  at  small  scales  must  be  available  if  these  scales  are  to  be  initialized  in  a 
circulation  model.  For  scales  comparable  to  or  larger  than  the  Rossby  radius,  description  of 
the  mass  field  is  sufficient  information,  and  the  velocities  adjust  to  reach  geostrophic 
equilibrium  with  the  observed  mass  field.  Selecting  an  initial  velocity  field  is  not  critical 
for  these  scales;  the  difference  between  starting  with  geostrophic  initial  velocities  or  zero 
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initial  velocities  is  reflected  only  in  some  damping  of  the  density  field  in  the  latter  case, 
because  more  energy  is  extracted  from  the  initial  disturbance  to  spin-up  the  velocity  field. 

The  DFI  results  reveal  meanders  with  horizontal  length  scales  of  20  to  25  km;  to 
determine  whether  they  are  influenced  by  the  specification  of  the  initial  velocity  field  or 
not,  an  estimate  of  the  Rossby  radius  is  needed.  In  a  continuously  stratified  fluid,  with  a 
vertical  scale  much  smaller  than  the  horizontal  scale  and  no  variations  in  bottom 
topography,  the  normal  modes  solution  to  the  motion  can  be  applied,  based  on  a  separation 
of  variables  between  depth  independent  horizontal  information  and  horizontally 
independent  vertical  structure.  Each  normal  mode  has  its  own  vertical  structure  and 
behaves  like  an  equivalent  homogeneous  system  with  depth  He  (equivalent  depth). 

The  first  mode  (or  mode  zero)  is  the  barotropic  mode  with  equivalent  depth  equal 
to  the  total  depth  of  the  water  column.  For  systems  predominantly  baroclinic,  the  vertical 
structure  of  the  higher  order  modes  is  more  useful.  The  normal  modes  (barotropic  and 
baroclinic)  can  be  found  by  solving  the  following  eigenvalue  problem  (Kundu,  1990): 


dz \N2dz  j 
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where  N2  is  the  depth  dependent  buoyancy  frequency,  is  the  vertical  structure  function 
of  any  variable,  cn2  is  a  separation  constant  with  units  of  velocity  (reciprocal  of  the 
eigenvalues)  and  n  is  the  mode  number.  This  equation  was  solved  for  the  normal  modes  of 
vertical  velocity,  using  an  average  profile  for  the  buoyancy  frequency  computed  from  all 
the  CTD  casts  in  the  cruise  assuming  a  maximum  common  depth  of  2000  m,  with 
specification  of  appropriate  boundary  conditions  for  the  vertical  function  (rigid  lid  at  the 
surface  and  bottom): 

yn  =  0,  z  =  0,—H.  (4.7) 

The  modal  shapes  for  the  first  four  baroclinic  modes  for  vertical  velocity  are  shown 
in  Figure  14.  The  first  mode  contains  about  20%  of  the  total  fit  distributed  by  all  the 
resolved  modes  (250).  The  other  modes  contribute  less  than  6%  to  the  total  fit.  A 
comparison  of  Figure  14  with  Figure  1 1  shows  that  the  first  mode  describes  well  the  region 
of  maximum  vertical  velocities  between  50  and  300  m;  a  combination  of  the  remaining 
modes  can  account  for  the  maximum  extension  of  the  eddy,  varying  from  500  to  about 
1000  m. 
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Figure  14.  Normalized  modal  shapes  for  vertical  velocity  computed  from  the  mean 
stratification. 


Each  of  the  resolved  modes  has  a  corresponding  eigenvalue  (k=l/cn2)  from  which 
the  associated  internal  Rossby  radius  can  be  computed: 


Table  3  lists  the  results  for  the  first  three  baroclinic  modes. 


(4.8) 


Mode  No. 

X 

Cn  (m/s) 

Rd(km) 

1 

0.235 

2.06 

23 

2 

1.070 

0.97 

11 

3 

2.474 

0.64 

7 

Table  3.  Characteristics  of  the  First  Three  Baroclinic  Modes 


The  first  internal  or  baroclinic  Rossby  radius  of  deformation  is  the  horizontal  length 
scale  associated  with  coastal  phenomena,  such  as  upwelling  and  meandering  jets,  and 
eddies  (Gill,  1982).  These  results  are  in  good  agreement  with  the  estimates  of  the  Rossby 
radius  off  the  coast  of  Oregon,  given  in  Kundu  et  al.  (1975),  obtained  with  a  similar  method 
to  the  one  described  above  (decomposition  into  normal  modes). 

It  follows  from  the  above  that  it  is  reasonable  to  assume  that  the  vertical  motion 
associated  with  the  meanders  observed  in  the  mean  flow  has  a  vertical  structure  dominated 
by  the  first  mode,  which  has  a  Rossby  radius  of  23  km.  The  length  scales  of  the  up-  and 
downwelling  patterns  shown  in  Figure  9  are  comparable  to  the  Rossby  radius.  Thus,  during 
the  model  initialization,  it  is  the  velocity  field  that  primarily  adjusts  to  the  observed  density 
field  at  the  scales  of  the  meanders. 

Therefore,  the  choice  of  the  initial  condition  for  the  velocity  field  is  not  critical  to 
achieve  the  same  qualitative  results  from  the  model  (meandering  of  the  mean  flow,  with 
associated  patches  of  rising  and  sinking  motion).  This  is  illustrated  in  Figure  15,  which 
shows  the  vertical  velocity  at  100  m  depth,  obtained  with  zero  vertical  shear  in  the  initial 
velocity  field  as  opposed  to  the  geostrophic  shear  used  to  obtain  the  results  shown  earlier. 
Comparing  this  figure  with  Figure  9,  it  is  seen  that  the  same  patterns  of  vertical  motion  are 
found  with  the  different  initial  velocity  field,  as  expected.  Another  result  that  is  consistent 
with  the  considerations  above  is  the  weaker  magnitudes  of  the  motion,  which  are  decreased 
by  more  than  50  %  at  this  level  relative  to  the  results  in  Figure  9.  Figures  16  and  17  show 
that  the  vertical  extent  of  the  vertical  motion  cells  is  also  somewhat  smaller  when  the  initial 
velocity  field  has  no  shear.  The  same  qualitative  agreement  with  the  corresponding  Figures 
1 1  and  12  is  clear,  but  the  regions  of  rising  and  sinking  motion  now  show  weaker  maximum 
and  minimum  vertical  velocities  in  each  section. 

The  question  then  arises  as  to  which  set  of  initial  conditions  yields  results  closer  to 
reality,  from  a  quantitative  point  of  view.  The  results  obtained  when  the  initial  velocity  field 
is  assumed  geostrophic  show  a  better  agreement  with  observed  currents.  To  set  the  initial 
shear  velocities  to  zero  is  equivalent  to  saying  that  no  a  priori  knowledge  is  available  about 
the  expected  velocities  in  the  eddy.  However,  scaling  arguments  can  show  that  the  eddy 
should  be  in  a  quasigeostrophic  balance  with  the  observed  mass  field  (the  Rossby  number 
is  of  order  0.1),  therefore  since  no  detailed  observations  of  the  initial  velocity  field  are 
available,  assuming  initial  geostrophically  balanced  currents  seems  the  better  assumption. 
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(alongshore)  (Min, Max)  W  =  (-6.61,7.304)  m/day 


3.  Data  related  problems 


The  underlying  assumption  in  the  data  analysis  described  in  section  ELB  is  that  the 
CTD  stations  (plus  the  four  SAIL  stations  and  the  satellite  data)  are  synoptic.  This  is  a 
common  assumption  in  oceanography,  perhaps  the  only  viable  solution  to  the  financial 
constraints  and  the  difficulty  of  data  collection  in  the  ocean.  Other  authors  who  have 
studied  meso-  and  submesoscale  phenomena  (Fiekas  et  al.,  1994;  Strass,  1994;  Pollard  and 
Regier,  1992;  Le  Groupe  Tourbillon,  1983)  make  no  particular  point  of  the  fact  that  this 
basic  assumption  is  made,  even  though  the  observations  are  collected  over  a  number  of 
days  or  even  a  week. 

The  scales  of  the  phenomena  being  studied  are  related  to  the  more  or  less  validity 
of  the  synopticity  assumption.  Data  collected  for  large  scale  phenomena  with  Rossby 
numbers  much  smaller  than  unity  or  time  scales  much  larger  than  the  duration  of  an 
hydrographic  cruise  can  be  considered  synoptic  without  any  major  argument  to  the 
contrary.  When  the  scales  being  studied  are  comparable  to  the  duration  of  the  cruise  or 
characteristic  length  scales  of  other  signals,  then  the  signal-to-noise  ratio  in  the 
measurements  becomes  smaller,  and  aliasing  and/or  deformation  of  the  scales  under 
analysis  may  occur. 

The  CTD  data  in  this  cruise  took  about  3  days  to  collect.  Figure  18  shows  the 
cumulative  time  interval  between  the  beginning  of  each  cast.  The  difference  in  time 
between  any  two  stations  is  found  by  subtracting  the  numbers  between  these  stations  and 
the  immediately  preceding  stations  (for  instance,  the  time  interval  between  stations  16  and 
15  is  40.9  -  21.6  -  19.3  hours.  Note  that  the  CTD  number  sequence  -Figure  2-  does  not 
correspond  to  the  order  in  which  the  stations  were  occupied).  The  worst  case  of  possible 
aliasing  or  distortion  occurs  when  stations  are  close  in  space  but  separated  by  several  hours 
in  time.  For  example,  stations  16  and  15  are  apart  by  about  20  km  in  space  and  19  hours  in 
time.  Thus  an  oscillation  in  time,  uniform  over  the  entire  grid,  with  a  period  of  38  hours  has 
the  potential  to  explain  all  of  the  differences  between  the  observed  values  at  stations  15  and 
16.  Note  that  the  spatial  aliasing  problem  in  this  case  is  solved  by  sampling  theory,  i.e., 
scales  smaller  than  the  separation  between  stations  are  not  resolved. 

This  section  examines  two  possible  sources  of  error  in  the  data  used  in  this  study 
and  tries  to  draw  some,  if  not  quantitative,  at  least  qualitative  conclusions  on  their  effect  on 
the  results  presented  in  the  preceding  sections.  These  possible  sources  of  error  are  a  direct 
consequence  of  the  lack  of  exact  synopticity  in  the  hydrographic  data. 
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Distance  Alongshore  (km) 


Figure  18.  Cumulative  time  interval  between  CTD  stations,  in  hours.  Refer  to 
Figure  2  for  stations  numbers  {not  related  to  the  order  the  stations  were  occupied). 

a.  Internal  tides 

The  interaction  of  the  barotropic  tidal  motions  with  the  coastline  at  the  shelf 
or  slope  causes  a  vertical  oscillation  of  the  pycnocline  with  the  same  range  of  frequencies 
of  the  tides.  These  oscillations  propagate  away  from  the  region  of  generation  as  internal 
waves  (internal  tides);  the  amplitude  of  the  internal  tides  depends  on  the  strength  of  the 
forcing,  which  in  turn  is  a  function  of  the  bottom  slope,  the  barotropic  tide-induced  across 
shelf  velocity  and  the  stratification  of  the  water  column  (Baines,  1986).  Their  phase  speed 
is  a  function  of  the  buoyancy  (TV)  and  inertial  if)  frequencies. 

Internal  waves  can  only  exist  in  the  range  of  frequencies 

f<U)<N.  (4.9) 
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Hence,  the  main  components  of  the  barotropic  tide  that  are  forcing  internal 
tides  are  semidiurnal.  Several  authors  (e.g.,  Rosenfeld,  1990)  have  found  that  the  frequency 
associated  with  the  largest  variance  in  internal  oscillations  close  to  the  shore  is  the  one  that 
corresponds  to  the  M2  component  (about  0.08  cycles/hour). 

Two  theories  try  to  explain  the  propagation  of  internal  waves.  One  approach 
uses  ray-like  propagation,  and  explains  important  phenomena  such  as  reflections  from  the 
bottom.  Depending  on  a  critical  slope  angle  determined  by  the  wave  characteristic 
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the  wave  will  be  either  transmitted  or  reflected  (Pond  and  Pickard,  1983).  This  can 
account  for  the  concentration  of  energy  at  some  depths,  as  waves  are  forced  to  remain  at 
selected  levels.  Another  approach  decomposes  the  oscillation  into  an  infinity  of  vertical 
dynamic  modes  similar  to  the  ones  depicted  in  Figure  14,  and  can  explain  the  difference  in 
the  amplitudes  of  the  oscillations  at  different  depths.  Usually  the  first  baroclinic  mode 
describes  most  of  the  variance  (Rosenfeld,  1990),  but  other  researchers  have  found  that  it 
can  take  a  combination  of  several  of  the  lowest  modes  to  explain  the  vertical  structure  of 
the  internal  tide  (Chiswell,  1994). 

In  any  case,  internal  wave  dynamics  are  dependent  on  such  a  variety  of 
contributions  that  the  result  is  a  fairly  complex  system  that  may  or  may  not  be  causing 
problems  in  the  CTD  data  sampling  in  this  study.  The  motivation  to  look  at  this  problem 
came  from  the  finding  of  appreciable  vertical  oscillations  in  the  upper  500  m  of  the  water 
column,  based  on  CTD  casts  taken  during  an  anchor  station  about  350  km  off  Pt.  Arena, 
thus  much  farther  offshore  than  the  present  study  region,  in  August  1993  (Michael  Kosro, 
personal  communication).  These  oscillations  had  a  well  marked,  vertically  coherent 
semidiurnal  character,  with  vertical  displacements  of  the  thermocline  of  about  10-15  m, 
(30  m  oscillations  in  the  isotherms  at  300  m  depth).  Such  a  large  and  organized  oscillation 
at  the  semidiurnal  period,  if  present  in  the  study  region,  would  undoubtedly  affect  the 
observations  in  the  CTD  data  used  in  this  study. 

The  only  time  series  available  in  this  study  are  from  the  current  meter 
records.  The  instruments  recorded  velocity  and  temperature  at  several  depths  (section 
H.A.2).  Unfortunately,  the  vertical  resolution  is  very  poor,  and  the  period  covered  does  not 
include  the  time  of  the  CTD  sampling.  In  the  following  analysis,  therefore,  only  the  last 
week  of  current  meter  data  were  used  (spanning  from  12  to  20  May  1993).  Figure  19  shows 
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Figure  19.  Predicted  sea  level  at  the  San  Francisco  bay  entrance,  for  May  1993,  in 
meter.  The  dash-dot  vertical  lines  mark  the  period  of  the  CTD  observations. 


the  sea  level  change  as  predicted  from  tide  tables  for  the  San  Francisco  Bay  entrance  for 
May  1993;  it  can  be  seen  that  the  barotropic  M2  forcing  is  roughly  the  same  during  the  CTD 
cruise  period  and  the  week  before  (the  CTD  period  was  coincident  with  a  new  moon  phase, 
thus  with  small  M2  amplitudes). 

A  classical  least-squares  harmonic  analysis  was  carried  out  on  each  of  the 
18  current  meter  time  series  to  evaluate  the  spatial  phase  and  amplitude  distribution  of  the 
most  important  tidal  constituents;  by  doing  this  at  all  levels  and  comparing  the  phase  and 
amplitude  differences,  an  idea  of  the  modal  distribution  can  be  obtained.  The  harmonic 
analysis  takes  the  time  series  and  fits  amplitude  and  phase  to  sinusoidal  basis  functions  with 
known  frequencies.  With  these  short  time  series,  only  a  few  components  are  resolved. 
Longer  time  series  are  not  advised  because  of  the  relatively  short  temporal  scale  of  the 
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internal  tides.  In  most  of  the  cases  (1 1  out  of  18  current  meters)  the  largest  amplitude  was 
associated  with  the  M2.  In  the  remaining  cases,  the  largest  amplitude  was  either  located  at 
wave-wave  components,  or  at  frequencies  that  are  not  supported  by  internal  waves, 
without  any  dominant  component. 

Figure  20  shows  the  upper  500  m  of  an  alongshore  vertical  cross  section  of 


Figure  20.  Alongshore  vertical  cross  section  of  temperature  (°C)  from  raw  CTD 
data.  The  horizontal  axis  is  the  station  number  in  the  section. 


raw  temperature  data  (not  interpolated).  In  the  vicinity  of  the  eddy,  the  temperature 
differences  observed  between  CTD  stations  are  (XIO*1  °C),  being  larger  in  the  lower  levels. 
The  fitted  amplitudes  for  the  M2  were  in  general  one  order  of  magnitude  smaller. 

The  attempt  to  fit  normal  modes  to  the  M2  amplitudes  was  inconclusive  due 
to  the  poor  vertical  resolution  of  the  continuous  records.  The  horizontal  distribution  of 
fitted  phases  appeared  random,  without  any  distinct  pattern.  The  assumption  of  random, 
rapidly  varying  phase  distribution  of  the  baroclinic  tides  has  already  been  used  by 
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Noble  et  al.  (1987)  in  their  analysis  of  tidal  currents.  The  results  above  seem  to  support  this 
assumption  in  the  present  situation. 

The  cross-correlation  between  temperature  and  horizontal  velocity  time 
series  at  each  instrument  was  also  carried  out  to  try  to  identify  internal  wave  characteristics 
at  the  M2  frequency,  using  cross-spectra  techniques  described  by  many  data  analysis 
references  (e.g.,  Bendat  and  Piersol,  1993).  According  to  Figure  21,  the  amplitude  of  the 
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Figure  21.  Relationship  between  the  phase  and  group  velocities  and  the  velocity 
vector  in  a  propagating  internal  wave  (the  y-axis  is  into  the  paper).  C  is  the  phase 
velocity,  Cg  is  the  group  velocity,  k  and  m  are  the  wavenumbers  (from  Kundu, 
1990). 

temperature  disturbance  associated  with  internal  waves  is  90°  out  of  phase  with  the 
u-component  of  velocity  and  in  phase  with  the  v-component.  The  results  of  the  analysis  of 
the  coherence  function  and  phase  spectra  for  all  the  instruments  are  listed  in  Table  4.  An 


Station  num¬ 
ber  (depth) 

Tvs.  U 

Tvs.  V 

Coherence 

Phase 

Coherence 

Phase 

13  (100) 

0.90 

95 

0.90 

175 

9(100) 

0.80 

160 

0.80 

85 

7(100) 

0.70 

170 

0.85 

95 

8  (100) 

0.90 

130 

0.90 

150 

3 (150) 

0.50 

135 

0.4 

55 

7 (150) 

0.05 

70 

0.15 

45 

Table  4.  Results  of  Currentmeter  Cross-correlation  Analysis  for  the  M2 
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Table  4.  Results  of  Currentmeter  Cross-correlation  Analysis  for  the  M2  (Continued) 

inspection  of  the  results  shows  that  in  only  a  couple  of  instruments  are  the  temperature  and 
velocities  well  correlated  at  the  M2  frequency,  with  phase  values  that  match  the 
requirements  above  (station  13  at  100  m  and  station  3  at  300  m,  with  an  allowance  made 
for  phase  uncertainty).  A  variety  of  situations  occur  for  other  instruments:  relatively  high 
values  of  coherence  are  associated  with  phase  values  that  seem  reversed  (e.g.,  stations  9  and 
7  at  100  m),  or  relatively  high  values  of  coherence  are  associated  with  meaningless  phase 
values  (e.g.,  station  13  at  150  m,  station  8  at  600  m),  or  temperature  and  velocities  are 
hardly  or  not  correlated  at  all  at  this  frequency  (e.g.,  station  7  at  150  m,  station  9  at  300  m). 
Moreover,  the  correlation  and  phase  values  are  not  consistent  over  an  appreciable  depth  of 
the  water  column.  These  results  show  a  complicated  vertical  and  horizontal  pattern  of 
waves  at  the  M2  frequency  whose  influence  in  the  sampling  is  very  likely  random. 

The  vertical  correlation  along  the  same  water  column  was  further 
investigated  by  computing  the  correlation  and  coherence  functions  between  the 
temperature  time  series  from  the  instruments  in  the  same  mooring.  The  values  found  at  all 
the  moorings  were  very  small,  even  between  instruments  separated  by  only  50  m.  The 
coherence  values  at  the  M2  frequency  were  generally  smaller  than  0.5.  The  vertical 
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extension  of  the  vertical  motion  patterns  computed  by  the  model  and  presented  in  section 
IV.  A.  1  cannot  therefore  be  explained  by  an  eventual  aliasing  of  internal  tides  in  the  CTD 
data.  Internal  tides  may  be  present,  as  suggested  by  the  results  of  the  harmonic  and  cross¬ 
correlation  analyses,  but  the  spatial  distribution  of  phases  and  amplitudes  seems  random 
enough  to  cause  small  contributions  to  the  field  captured  in  the  hydrographic  data.  The 
chances  that  these  randomly  varying  fields  would  become  coherent  and  influence  the 
sampling  in  the  worst  possible  manner  seem  very  small. 

Several  factors  may  be  contributing  to  this  random  distribution.  The  period 
of  the  CTD  sampling  was  during  a  time  of  weak  amplitudes  in  the  barotropic  semidiurnal 
frequency  range,  thus  the  forcing  is  minimal.  Between  22  and  25  May  there  was  a  mix  of 
relatively  strong  wind  events  and  relaxations  (Figure  3),  causing  the  stratification  to  vary 
and  affecting  the  maintenance  of  internal  oscillations.  Finally,  the  area  where  this  study 
took  place  may  be  sufficiently  close  to  the  generation  sites  of  the  internal  tides  (composed 
of  an  infinity  of  source  points  along  an  irregular  shelf  or  slope)  for  the  internal  wave  field 
to  be  spatially  incoherent  and  appear  random.  Further  offshore  the  wave  field  would 
organize  itself  into  a  vertically  and  horizontally  coherent  set  of  modes  of  propagation, 
which  could  explain  the  measurements  of  M.  Kosro  away  from  the  coast. 

Some  questions  concerning  the  propagation  of  internal  tides  were  raised  by 
a  study  of  well  defined  internal  tide  oscillations  700  km  west  of  Santa  Barbara,  California 
(Levine  and  Richman,  1989).  By  performing  an  elliptical  decomposition  of  the  observed 
signals,  they  found  that  the  predominant  internal  waves  in  that  area  over  a  period  of  a  few 
days  propagated  in  a  direction  parallel  to  the  coast,  contradicting  the  theories  that  predict 
propagation  normal  to  the  coast  (Baines,  1986).  They  suggest  that  internal  tides  can  be 
refracted  by  the  density  structure  of  the  ocean  and  thus  propagate  in  a  variety  of  directions. 
Since  the  density  structure  in  the  eddy  in  the  present  study  is  fairly  complex,  these  findings 
support  the  present  assumption  of  small  spatial  coherence  and  highly  random  nature  of  the 
internal  tide  signal  with  a  minimal  influence  on  the  observations. 

b.  Eddy  rotation 

This  subsection  is  devoted  to  an  analysis  of  the  effect  of  eddy  rotation  on  the 
solutions  shown  earlier.  This  is  perhaps  a  more  acute  problem  than  the  internal  tides 
problem,  because  there  is  no  way  to  quantify  the  effect  of  sampling  a  rotating  body  of  water 
while  the  data  is  being  collected.  The  analysis  presented  below  assumes  that  there  is  a 
wavenumber  in  the  eddy  at  least  equal  to  two;  it  will  be  shown  that  a  more  reasonable 
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estimate  is  at  least  four,  based  solely  on  the  time  intervals  between  stations.  The  same 
reasoning  is  then  applied  to  show  that  the  DFI  solution  of  a  wavenumber  five  (Figure  7)  is 
very  reasonable.  Translation  of  the  eddy  is  neglected.  According  to  Nof  (1981),  considering 
a  well  mixed  core  300  m  deep,  the  translation  speed  would  be  about  1.4x1  O'2  m/s  westward, 
giving  only  3  to  4  km  of  offshore  displacement  during  the  three  days  of  CTD  data 
collection. 

The  rotation  of  the  eddy  can  be  estimated  by  taking  40  cm/s  as  being 
representative  of  the  tangential  horizontal  speed  at  30  to  40  km  from  the  center  at  the 
surface,  and  assuming  a  solid  body  rotation  for  simplification.  Under  these  premisses  the 
eddy  takes  about  6  days,  or  twice  the  period  of  the  CTD  sampling,  to  complete  one  rotation. 
The  problem  here  is  to  determine  how  the  sub-eddy  scale  features  have  been  sampled 
during  the  half  eddy  rotation.  Clearly,  the  features  themselves  will  not  propagate  at  the 
same  speed  of  the  almost  circular  mean  flow,  because  there  exists  a  potential  vorticity 
gradient  across  the  mean  flow  that  acts  like  an  equivalent  beta  effect  for  Rossby  waves,  thus 
retarding  the  propagation  speed  of  the  features  relative  to  the  mean  flow. 

A  theoretical  estimate  of  the  phase  speed  of  the  sub-eddy  meanders  can  be 
obtained  by  applying  the  results  of  baroclinic  instability  theory  in  the  special  case  of  a 
barotropic  basic  state  (Holton,  1979,  pp.  216),  with  the  beta  effect  being  replaced  by  its 
equivalent  potential  vorticity  gradient,  which  amounts  to  the  gradient  of  relative  vorticity 
across  the  mean  flow  due  to  the  flat  bottom  assumption.  The  phase  speed  of  the 
disturbances  is  therefore  estimated  from: 

c-u„- (^z*2)  «■”> 


where  Um  is  the  mean  speed  of  the  layer,  k  is  the  azimuthal  wavenumber,  n  is  the  radial 
coordinate  in  the  natural  coordinate  system  and  C,  is  the  azimuthally  averaged  relative 
vorticity  given  by  (in  natural  coordinates): 


r  =  Y.-W 
**  R  9n 


(4.12) 


R  being  the  distance  to  the  center  of  curvature  and  V  the  tangential  speed.  The  effect  of  the 
second  term  on  the  right  hand  side  of  Equation  4.11  is  expected  to  be  small  for  high 
wavenumbers  (phase  speeds  of  about  90  to  95%  of  the  mean  azimuthal  flow),  in 
accordance  with  the  results  from  the  numerical  experiments  of  McCalpin  (1987). 
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The  mean  tangential  velocities  were  computed  as  a  function  of  distance  to 
the  center  of  the  eddy,  by  taking  azimuthal  averages  of  horizontal  speed  at  several  radial 
distances  from  the  center,  and  then  averaging  in  the  vertical  over  the  first  500  m,  i.e.,  over 
the  mean  vertical  extent  of  the  sub-eddy  features.  The  results  show  mean  tangential  speeds 
of  about  6  cm/s  at  10  km  from  the  center  of  the  eddy,  increasing  to  about  14  cm/s  at  30  km 
and  17  cm/s  at  35  km,  and  then  decreasing  to  about  12  cm/s  at  50  km  from  the  center.  Using 
these  results  in  Equations  4.1 1  and  4.12,  with  some  radial  averaging  over  the  radial  scale 
of  the  disturbances,  gives  phase  speeds  for  the  sub-eddy  features  of  about  14  cm/s  at  a 
distance  of  30  km  from  the  center  of  the  eddy,  that  is,  about  16  days  to  complete  one 
rotation. 

This  estimate  for  the  movement  of  the  meanders  agrees  well  with  that 
observed  during  the  model  integrations;  the  animation  of  the  horizontal  velocity  field  at 
each  time  step  showed  an  angular  displacement  of  the  sub-eddy  features  of  about  10°  about 
the  center  in  12  hours,  or  18  days  for  a  complete  rotation. 

Therefore,  during  the  3  days  of  CTD  data  collection,  the  meanders  are 
estimated  to  have  moved  about  60°  cyclonically.  It  is  clear  that  a  wavenumber  two  around 
the  eddy  cannot  create  the  appearance  of  a  wavenumber  five  seen  in  the  results  of  section 
IV. A.  1  by  aliasing  only,  because  the  same  trough  (or  ridge)  could  not  be  sampled  three 
times  during  the  cruise.  Similarly,  an  analysis  of  Figure  18  shows  that  there  is  no  time 
interval  between  CTD  stations  large  enough  to  cause  aliasing  of  wavenumbers  smaller  than 
the  wavenumber  five  found  in  the  DFI  solution.  Wavenumbers  larger  than  five  can 
theoretically  be  aliased,  but  this  would  imply  that  the  length  scale  of  the  sub-eddy  features 
would  have  to  be  smaller  than  the  Rossby  radius,  which  is  the  preferred  length  scale  for  the 
growth  of  baroclinic  instabilities  in  a  /-plane  (Gill,  1982).  Thus  aliasing  in  the  data  due  to 
meander  propagation  is  not  a  cause  for  errors  in  the  solution. 

Comparing  the  time  interval  between  CTD  stations  in  Figure  18  with  the 
streamline  and  vertical  motion  patterns  in  Figure  9,  and  using  the  new  information  on  the 
propagation  of  these  sub-eddy  features,  it  is  seen  that  the  maximum  time  interval  between 
adjacent  stations  is  34  hours  (stations  20  and  21);  these  two  stations  are  located  at  the  edges 
of  a  sinking  cell  that  lays  upstream  of  a  relatively  broad  region  of  upwelling.  Hence,  station 
20  sampled  the  leading  edge  of  the  cyclonically  rotating  downwelling  cell  and  station  21 
sampled  the  trailing  edge  of  the  same  cell,  34  hours  later  when  the  cell  had  already  rotated 
about  30°.  Therefore,  the  net  effect  of  the  asynoptic  sampling  in  this  instance  is  to  decrease 
the  apparent  width  of  the  downwelling  cell.  The  same  reasoning  applies  for  stations  19  and 
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22  and  others.  In  a  global  sense,  this  explains  the  decrease  in  the  horizontal  scale  of  the 
features  sampled  against  the  sense  of  rotation  of  the  eddy,  as  exemplified,  and  the  increase 
in  the  scales  of  the  features  sampled  along  the  sense  of  rotation,  such  as  the  sequences  of 
stations  28  to  21. 

Other  sequences  of  stations  are  not  significantly  affected  by  the  propagation 
of  the  meanders  because  the  time  between  stations  was  not  long  enough  to  cause  an 
appreciable  change  in  length  scales;  these  include  the  sequence  of  stations  8,  7,  6, 11,  12 
and  13  which,  although  sampled  along  the  rotation,  are  apart  by  only  about  13  hours  total. 

It  then  seems  reasonable  to  conclude  that  the  effect  of  the  eddy  rotation,  or 
meander  propagation,  in  the  CTD  data  is  to  distort  the  horizontal  scales  of  some  of  the 
features,  namely  those  located  on  the  southeast  side  of  the  eddy  (right  side  in  Figure  7).  The 
features  on  the  northwest  side  (left  side  in  Figure  7)  must  have  maintained  their  scales  with 
minimum  distortion  because  the  time  interval  between  stations  is  relatively  short  compared 
to  the  angular  period  of  the  meanders.  Since  the  troughs  and  ridges  associated  with  the  sub¬ 
eddy  scales  are  about  70°  apart,  it  would  require  a  sequence  of  stations  with  time  intervals 
on  the  order  of  several  days  to  cause  aliasing  of  this  wavenumber.  Thus,  although 
wavenumber  five  is  deformed,  it  is  an  alias-free  estimate  of  the  measured  field. 

c.  Interpolation  parameters 

The  data  analysis  described  in  section  n.B  makes  use  of  the  MQ-B 
interpolation  method  to  fill  the  model  domain  with  the  CTD  and  satellite  data.  This  method 
makes  use  of  several  interpolation  parameters,  the  most  important  of  which  are  the 
interpolation  constant  c  and  the  exaggeration  applied  to  the  vertical  axis  (see  section 
D.B.l).  In  addition,  it  was  also  described  in  section  n.B  that  the  analysis  procedure  used  for 
the  temperature  and  salinity  fields  was  not  the  same,  because  the  number  of  data  points  for 
temperature  and  salinity  are  different. 

At  that  point,  the  criterion  for  selecting  the  interpolation  parameters  was  to 
achieve  the  best  fit  to  the  observed  profiles  (checked  by  plotting  each  observed  profile 
against  several  interpolated  fields)  and  to  assure  a  minimum  difference  between  extreme 
data  values  and  extreme  interpolated  values.  As  noted,  the  application  of  these  criteria 
along  with  the  constraint  in  the  number  of  observations  resulted  in  the  choice  of  different 
sets  of  parameters  for  temperature  and  salinity. 

The  process  of  determining  the  sensitivity  of  the  circulation  model/DFI 
results  to  the  set  of  parameters  chosen  requires  two  steps:  the  interpolation  of  the  data  with 
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different  parameters,  yielding  slightly  different  analyzed  fields  each  time;  and  running  the 
model  with  each  of  the  newly  analyzed  fields  to  compare  results.  Several  such  runs  were 
made.  It  was  noticed  that  although  the  analyzed  fields  were  very  sensitive  to  the 
interpolation  constant  c,  the  dynamically  balanced  velocity  fields  maintained  the  same 
qualitative  agreement  in  the  location  of  the  meanders  and  vertical  motion  areas;  the 
differences  were  found  in  the  magnitude  of  the  dynamically  balanced  fields  (either  stronger 
or  weaker  velocities). 

Another  sensitivity  test  was  made  to  estimate  the  influence  of  analyzing 
temperatures  and  salinities  differently.  In  this  case,  the  two-dimensional  interpolation  of 
the  CTD  data  at  each  level  (see  section  n.B.2.a)  was  also  performed  for  the  salinities,  and 
the  result  was  then  three-dimensionally  interpolated  into  the  model  domain  using  the  same 
parameters  as  for  the  temperature  interpolation.  The  result  was  a  salinity  field  that  looked 
better  but  had  worse  fits  with  the  observed  profiles  and  maxima  and  minima.  The  model 
DF1  solution  showed  a  weaker  field  of  vertical  motion,  but  in  good  qualitative  agreement 
with  the  original  results. 

These  considerations  demonstrate  the  more  or  less  dependence  of  the  results 
on  the  MQ-B  method  because  of  the  limited  data  available.  The  choice  of  the  parameters 
utilized  in  this  study  was  dictated  by  criteria  that  guaranteed  (after  a  long  process  of  trial 
and  error)  the  closest  possible  match  between  observed  and  analyzed  fields  in  the  regions 
where  data  was  available.  Since  the  computed  velocities  and  the  observed  currents  from  the 
moorings  agree  much  better  for  the  fields  analyzed  as  described  in  II.B  than  they  do  for  the 
other  analysis  made  for  sensitivity  purposes  (again  with  a  limited  set  of  data  for 
comparison),  it  is  felt  that  the  analyzed  fields  from  n.B  and  the  associated  results  are 
correct  within  the  given  constraints. 

B,  QUASIGEOSTROPHIC  OMEGA  EQUATION  RESULTS 

The  vertical  motion  in  the  study  area  was  also  diagnosed  by  solving  the  Q-vector 
form  of  the  QG  omega  equation,  reproduced  below  for  convenience: 

2 

(N2V2+f0-^)w  =  2V»Q,  (4.13) 

Bz 

where 


60 
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with  zero  boundary  conditions,  as  described  in  section  ffl.B. 

The  right  hand  side  of  Equation  4.13  (known  as  the  forcing)  was  computed  from  the 
analyzed  density  field,  and  also  from  the  dynamically  balanced  density  field,  for 
comparison  purposes.  The  resulting  Q-vector  divergence  fields  were  smoothed  with  one 
pass  of  a  smoother-desmoother  filter  (Kreitzberg  and  Perkey,  1977;  Schuman,  1957),  as  is 
common  practice  in  similar  applications  (e.g.,  Pauley  and  Nieman,  1992).  An  example  of 
the  Q-vector  field  in  a  small  region  of  the  domain  at  100  m  depth,  obtained  from  the 
analyzed  density  field,  is  shown  in  Figure  22.  The  figure  shows  the  spatially  averaged 
Q-vectors,  dynamic  height  and  density  fields  in  the  area  of  the  meander  located  closest  to 
the  offshore  boundary  (lower  part  of  Figure  9;  this  figure  can  be  used  for  exact  location  of 
the  area  depicted  in  Figure  22).  Although  a  complete  analysis  of  the  dynamics  of  this 
meander  would  require  information  on  the  same  fields  above  and  below  the  selected  depth, 
Figure  22  provides  enough  information  to  assess  the  behavior  of  the  forcing  of  the  QG 
omega  equation  in  this  area.  A  trough  in  the  dynamic  height  field  extends  southward  and 
offshore  (roughly  from  the  point  (15,14)  to  the  point  (21,7)).  The  geostrophic  advections 
take  place  parallel  to  the  dynamic  height  contours.  Thus  upstream  of  the  trough  axis  there 
is  a  region  of  relatively  weak  cold  water  advection  from  the  center  of  the  eddy  and  negative 
vorticity  advection  from  the  upstream  ridge,  whereas  downstream  of  the  trough  there  is  a 
region  of  stronger  warm  water  advection  from  outside  of  the  eddy  and  positive  vorticity 
advection  from  the  trough.  Cold  and  warm  advection  by  the  geostrophic  currents  are  related 
to  downward  and  upward  motions,  respectively  (second  term  on  the  right  hand  side  of 
Equation  3.19).  Since  the  currents,  and  associated  advections,  generally  increase  toward  the 
surface  (the  first  term  on  the  right  hand  side  of  Equation  3.19),  upwelling  downstream  and 
downwelling  upstream  of  the  trough  should  be  expected.  Positive  and  negative  vorticity 
advection  are  responsible  for  the  propagation  of  the  meander  (in  the  absence  of 
divergence),  and  their  vertical  structure  also  contributes  to  the  vertical  motion. 

The  expected  patterns  of  upwelling  and  downwelling  cells  downstream  and 
upstream  of  the  trough,  respectively,  are  well  represented  by  the  Q-vectors,  which 
generally  point  toward  regions  of  rising  motion.  On  the  left  of  the  trough  the  vectors  are  in 
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Figure  22.  Field  of  Q-vectors  computed  from  the  analyzed  density  field  in  part  of  the 
domain,  at  100  m.  Also  shown  are  dynamic  height  (solid)  and  density  (dashed) 
contours.  The  X  and  Y  coordinates  are  the  same  as  in  Figures  9  and  23. 

general  oriented  away  (diverge)  from  the  region  of  cold  water  advection,  where  the 
dynamic  height  and  density  contours  intercept  at  a  greater  angle,  thus  denoting  the  area  of 
sinking  motion  upstream  of  the  trough.  Downstream  of  the  trough,  the  Q-vectors  are 
oriented  with  maximum  convergence  near  the  area  where  there  is  a  local  maximum  in  warm 
water  advection  (dynamic  height  and  density  contours  intercepting  almost  at  right  angles) 
and  consequent  upwelling.  Thus,  although  the  effects  of  differential  vorticity  advection 
have  not  been  considered,  such  an  analysis  is  consistent  with  a  similar  analysis  based  only 
on  the  changes  of  relative  vorticity  following  the  motion  of  a  parcel  through  the  meanders 
(section  IV.A.l). 

The  solution  for  w  on  the  left  hand  side  of  Equation  4.13  (known  as  the  response ) 
was  computed  with  a  simultaneous  over-relaxation  method  (Haltiner  and  Williams,  1980). 
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Figure  23  shows  the  field  of  vertical  motion  at  100  m  depth,  obtained  using  the  analyzed 
density  field.  A  very  good  qualitative  agreement  with  the  DFI  solution  depicted  in  Figure  9 
is  apparent.  The  two  solutions  reproduce  the  same  pattern  of  radially  aligned  patches  of 
sinking  and  rising  motion,  with  remarkable  precision  as  to  the  location  of  the  vertical 
motions  and  their  relative  magnitudes.  The  patterns  of  rising  and  sinking  motion  deduced 
above  from  the  Q-vectors,  dynamic  height  and  density  fields  (Figure  22)  near  the  meander 
in  the  lower  part  of  Figure  23  are  clearly  represented  and  in  agreement  with  the  forcing. 

Quantitatively,  the  vertical  velocity  solutions  in  Figures  9  and  23  are  somewhat 
different.  The  maximum  QG  vertical  motion  is  about  30%  larger  than  the  DFI  solution  in 
both  rising  and  sinking  areas.  This  difference  can  be  due  to  two  causes:  the  DFI  procedure 
is  appreciably  affecting  the  initial  density  field,  causing  the  resulting  DFI  vertical  velocity 
field  to  be  considerably  different  than  that  diagnosed  from  the  initial  density  field,  or/and 
the  difference  is  due  to  the  QG  assumptions.  The  relative  importance  of  the  former  can  be 
estimated  by  comparing  the  QG  omega  equation  results  obtained  from  the  analyzed  and 
from  the  dynamically  adjusted  density  fields.  Figure  24  shows  the  latter  results,  which 
should  be  compared  to  those  in  Figure  23.  It  is  seen  that  the  two  QG  solutions  are  not 
significantly  different,  that  is,  the  QG  vertical  motions  diagnosed  from  the  analyzed  and 
from  the  DFI  density  fields  are  almost  the  same.  This  result  is  explained  by  the  discussion 
of  the  geostrophic  adjustment  in  section  IV.A.2,  where  it  was  shown  that  the  scales  being 
initialized  in  the  model  are  comparable  to,  or  larger  than,  the  Rossby  radius.  Therefore,  the 
density  field  is  hardly  affected  by  the  DFI,  so  the  QG  solutions  should  not  be  very  different. 
Thus,  the  difference  between  the  two  solutions  for  w  (PE  and  QG)  is  entirely  due  to  the  QG 
assumptions. 

The  effects  of  the  QG  assumptions  in  the  QG  solution  for  w  can  be  assessed  from 
the  so-called  generalized  omega  equation,  which  contains  the  higher  order  effects  that  have 
been  neglected  in  the  QG  form,  namely  relative  vorticity,  tilting,  stability,  ageostrophic 
advection  and  vertical  advection  effects  (still  considering  adiabatic  solutions). 
Comparisons  between  the  generalized  and  the  QG  omega  equation  have  been  made  by 
Pauley  and  Nieman  (1992).  They  found  that  the  QG  omega  equation  solution  is  generally 
too  large  in:  1)  regions  of  strong  stability;  2)  regions  of  large  relative  vorticity;  3)  regions 
of  sinking  motion  in  cyclonic  curvatures  (rising  motion  in  cyclonic  flows  can  either  be 
under-  or  overestimated).  A  complete  analysis  of  the  generalized  omega  equation  is 
beyond  the  scope  of  this  work,  but  some  of  the  ageostrophic  effects  can  be  readily 
identified. 
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Figure  23.  Horizontal  cross  section  of  the  quasigeostrophic  vertical  velocity  field  at 
100  m  depth,  in  m/day,  from  the  solution  of  the  omega  equation  with  the  analyzed 
density  field.  The  orientation  and  horizontal  scales  are  the  same  as  in  Figure  9.  Note 
that  the  velocity  scale  is  different  from  Figure  9. 


The  stability  profile  used  in  the  QG  solution  in  Figure  23  is  a  spatial  average  of  all 
the  CTD  casts.  The  maximum  stability  is  found  around  30-40  m  and  below  that,  it 
decreases  sharply.  The  eddy  has  a  well  mixed  core  down  to  300  m,  thus  stability  must  be 
relatively  weak  in  the  eddy  at  100  m  depth. 

The  mean  relative  vorticity  of  the  eddy  is  about  0.1  /.  Thus  neglecting  relative 
vorticity  versus  planetary  vorticity  in  the  divergence  term  of  the  vorticity  equation  is  by 
itself  introducing  a  10%  error  in  the  response  of  the  QG  omega  equation.  This  effect  tends 
to  be  asymmetric,  causing  the  QG  solution  in  general  to  overestimate  w  in  regions  of 
positive  relative  vorticity  and  to  underestimate  w  in  regions  of  negative  relative  vorticity. 
This  effect  is  undoubtedly  contributing  to  some  of  the  overestimation  of  w  by  the  QG 
solution  in  the  meander  trough  shown  in  Figure  22. 

The  QG  solution  for  the  vertical  motion  presented  in  Figure  23  seems  to  be  strong 
in  both  rising  and  sinking  motions.  Thus  the  ageostrophic  advection  of  vorticity  and  density 
that  otherwise  would  be  opposing  the  geostrophic  advections  in  a  cyclonic  flow  are  missing 
in  the  whole  domain,  and  not  only  in  the  areas  of  sinking  motion.  Figure  25  shows  the 
ageostrophic  velocity  field  \ag,  computed  from: 


V  -  V 

v  DFI  v  g 


where  \DFI  is  the  horizontal  velocity  vector  resulting  from  the  DFI  procedure,  and  Vs 
the  initial  geostrophic  velocity,  as  well  as  the  density  and  dynamic  height  fields  for  the 
same  region  and  depth  as  in  Figure  22.  The  larger  vectors  shown  in  Figure  25  represent 
ageostrophic  velocities  that  are  only  about  10%  of  the  geostrophic  velocities  in  the  same 
area.  It  is  seen  that  upstream  of  the  trough,  the  ageostrophic  velocity  vectors  generally 
oppose  the  geostrophic  flow  and  are  associated  with  (weak)  warm  water  advection  and 
positive  geostrophic  vorticity  advection,  both  inducing  rising  motion.  This  region  is  where 
a  maximum  of  sinking  motion  in  the  QG  solution  is  related  to  the  cold  water  geostrophic 
advection.  Thus,  neglecting  ageostrophic  advection  of  density  and  geostrophic  vorticity, 
which  is  an  order  of  the  Rossby  number  effect,  causes  the  sinking  motion  to  be 
overestimated  upstream  of  the  trough  by  a  factor  of  about  10%. 

The  geostrophic  advection  of  ageostrophic  vorticity  upstream  of  the  trough,  which 
is  also  an  order  of  the  Rossby  number  effect,  can  be  estimated  by  a  careful  inspection  of 
Figure  25.  It  can  be  seen  that  the  ageostrophic  vectors  tend  to  curl  positively  in  two  main 
regions,  one  further  upstream  of  the  trough  and  another  just  downstream  of  the  trough. 
Between  these  two  regions  the  ageostrophic  vorticity  is  generally  negative.  Thus,  the 
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Figure  25.  Field  of  ageostrophic  velocity  vectors  at  100  m  depth,  for  the  same  part  of 
the  domain  as  in  Figure  22.  Also  shown  are  dynamic  height  (solid)  and  density 
(dashed)  contours. 

geostrophic  advection  of  ageostrophic  vorticity  is  positive  upstream  of  the  trough. 
Therefore,  neglecting  this  term  in  the  QG  solution  also  causes  an  overestimate  of  sinking 
motion  upstream  of  the  trough  by  about  10%. 

Downstream  of  the  trough,  the  ageostrophic  velocity  vectors  are  smaller  and  the 
effect  of  the  ageostrophic  advection  is  not  as  evident  as  it  is  upstream  of  the  trough.  In  a 
region  immediately  downstream  of  the  trough  axis,  the  ageostrophic  advection  is  weakly 
opposing  the  geostrophic  warm  water  advection  and  would  therefore  attenuate  the  QG 
estimate  for  rising  motion  there.  Moreover,  the  geostrophic  flow  is  advecting  negative 
ageostrophic  vorticity  into  the  same  area,  with  a  similar  contribution  for  the  attenuation  of 
the  QG  rising  motion  estimate.  Including  the  neglect  of  relative  vorticity  in  the  QG 


67 


divergence  term,  these  effects  all  indicate  an  overestimate  of  the  rising  motion  in  and 
immediately  downstream  of  the  trough  in  the  QG  solution.  Further  downstream  the  weak 
ageostrophic  advections  would  probably  slightly  enhance  the  rising  motion  computed  from 
the  QG  omega  equation,  and  other  effects  that  are  also  order  of  the  Rossby  number  (Pauley 
and  Nieman,  1992)  may  be  playing  a  role.  This  mixed  effect  is  in  accordance  with  the 
results  of  Pauley  and  Nieman  (1992)  for  regions  of  rising  motion  in  cyclonic  curvatures. 

Another  reason  that  could  account  for  some  of  the  differences  between  the  two 
solutions  is  the  high  order  numeric  computations  involved  in  the  determination  of  the 
forcing  of  Equation  4.13,  which  makes  the  divergence  of  the  Q-vectors  numerically 
sensitive  to  small  disturbances  in  the  density  and  geostrophic  velocity  fields. 

Although  the  analysis  above  was  focused  on  a  region  that  shows  the  maximum 
discrepancy  between  PE  and  QG  results,  it  can  be  generalized  to  the  other  meanders  as 
well.  It  then  seems  reasonable  to  conclude  that,  in  the  absence  of  other  effects,  it  is  a 
combination  of  eliminating  the  relative  vorticity  from  the  divergence  term  in  the  vorticity 
equation,  and  neglecting  the  ageostrophic  advection  and  vorticity,  that  is  contributing  the 
most  to  the  quantitative  differences  between  the  QG  and  the  DFT  estimates  of  vertical 
velocity,  with  the  QG  values  being  overestimated. 

The  vertical  structure  of  the  QG  vertical  velocity  solution  is  shown  in  Figures  26 
and  27.  These  figures  should  be  compared  to  the  DFI  solution  in  Figures  11  and  12, 
respectively.  Again,  a  very  good  qualitative  agreement  is  found  between  the  two  solutions, 
with  the  QG  solution  presenting  larger  values  for  maximum  and  minimum  velocities,  but 
with  some  good  agreement  in  other  areas.  The  vertical  extent  and  correlation  of  the  sinking 
and  rising  cells  of  vertical  motion  are  also  captured  by  the  QG  solution,  but  there  is  some 
variability  in  the  deeper  levels  that  is  not  found  in  the  DFI  solution.  The  area  close  to  the 
northeastern  boundary  (right  side  in  Figure  27)  where  the  CTD  profiles  were  extended  in 
depth  below  the  shallow  topography  also  shows  some  unrealistic  values  of  vertical  motion. 


Figure  26.  Vertical  cross  section  of  the  quasigeostrophic  vertical  velocity  field  in  the 
alongshore  direction,  in  m/day,  with  depths  in  meters.  The  section  shown  is  the  same 
as  in  Figure  1 1  (the  velocity  scales  are  different).  Note  that  the  vertical  axis  is  not 
linear. 
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X  index  (alongshore)  (Min, Max)  W  =  (-18.55,17.13)  m/day 


Y  index  (across  shore)  (Min, Max)  W  =  (-16.28,25.76)  m/day 


V.  CONCLUSIONS 


The  main  objective  of  this  study  was  to  achieve,  for  the  first  time,  a  quantitative 
description  of  the  three-dimensional  circulation  in  a  coastal  eddy,  including  the  vertical 
velocity.  The  temperature  and  salinity  data  collected  during  one  of  the  EBC  cruises  was 
thus  first  analyzed  using  multiquadric-biharmonic  interpolation.  The  resulting  density  field 
was  then  used  in  a  primitive  equation  model,  to  obtain  three-dimensional  fields  of 
dynamically  balanced  currents  and  densities  through  a  digital  filter  initialization  (DFI) 
technique.  The  DFI  is  equivalent  to  a  nonlinear  normal  mode  initialization,  and  thus  the 
dynamically  balanced  fields  resulting  from  the  model  integrations  under  DFI  are 
representative  of  the  slow  mode  at  the  time  of  the  observations. 

The  results  from  this  procedure  reveal  a  cyclonic  eddy  circulation  (relative  vorticity 
in  the  order  of  0.1-0.15/)  about  90-100  km  in  diameter,  with  a  total  of  five  meanders  each 
having  a  mean  horizontal  length  scale  of  about  25  km,  present  in  the  flow.  The  maximum 
horizontal  velocities  in  the  eddy  are  about  40-50  cm/s  at  the  surface,  in  good  agreement 
with  the  currents  observed  by  current  meters  moored  in  the  area,  and  with  the  hull  mounted 
ADCP  currents  estimates. 

The  meanders  are  associated  with  regions  of  rising  and  sinking  motion,  which  are 
shown  to  be  consistent  with  potential  vorticity  considerations  applied  to  a  water  parcel 
following  the  motion,  i.e,  downwelling  upstream  and  upwelling  downstream  of  the  troughs. 
These  patterns  of  vertical  motion  appear  to  be  radially  aligned,  distributed  on  either  side  of 
each  trough  axis,  and  extending  coherently  in  the  vertical  for  about  500  m  or  more.  The 
maximum  magnitudes  of  the  vertical  velocity  (either  rising  or  sinking)  occur  between  50 
and  150  m  depth,  and  are  about  20-22  m/day. 

The  vertical  distribution  of  the  vertical  velocity  amplitudes  is  fairly  well  described 
by  the  first  dynamical  mode,  which  has  a  characteristic  Rossby  radius  of  23  km.  This  was 
the  basis  for  a  discussion  on  the  process  of  geostrophic  adjustment  and  on  the  choice  of  an 
appropriate  initial  current  field.  It  was  shown  that  the  scales  determining  the  main  patterns 
of  vertical  motion  are  comparable  to  the  Rossby  radius,  and  thus  it  is  the  observed 
(analyzed)  density  field  that  determines  the  final  state  of  the  currents  at  these  scales.  The 
initial  velocity  field  was  therefore  assumed  to  be  geostrophic  while  the  depth  averaged  flow 
was  computed  assuming  a  level  of  known  motion  (2  cm/s  alongshore  to  the  northwest  at 
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2000  m  depth)  based  on  current  meter  data.  The  exact  magnitude  of  this  bottom  velocity 
was  found  not  to  be  significant  in  determining  the  final  dynamically  balanced  solution. 

At  this  stage  the  eddy  presents  a  closed  circulation  within  30  to  35  km  from  the 
center,  with  nearly  zero  average  vertical  velocities  within  this  area  at  all  levels.  This  means 
that  the  eddy  is  no  longer  acquiring  relative  vorticity  by  advection  from  other  sources,  and 
consequently  it  is  no  longer  spinning  up,  and  thus  it  is  about  to  detach  from  the  origin  region 
near  the  coast  and  to  drift  westward  (or  southwestward,  if  the  effect  of  the  California 
Current  farther  offshore  is  taken  into  account).  However,  the  effects  of  the  poleward 
undercurrent  close  to  the  shore  could  not  be  resolved  in  this  study  because  of  the  way  the 
data  were  processed  in  the  northeastern  part  of  the  domain.  Therefore,  allowance  must  be 
made  for  possible  vorticity  exchanges  in  the  regions  of  the  eddy  that  are  more  than  35  km 
away  from  the  center. 

Some  aspects  related  to  the  lack  of  exact  synopticity  in  the  hydrographic  data  were 
also  examined.  The  influence  of  the  internal  tides  on  the  sampling  is  thought  to  be  minimal, 
since  a  classical  harmonic  analysis  and  a  cross-correlation  analysis  of  the  current  meter 
time  series  indicate  that  the  spatial  phase  distribution  of  the  internal  tides  was  sufficiently 
random.  The  influence  of  the  eddy  rotation,  or  meander  propagation,  during  the  period 
covered  by  the  CTD  sampling  was  also  analyzed.  It  was  shown  that  a  wavenumber  five 
disturbance  in  the  mean  flow  around  the  eddy  could  not  have  been  aliased  during  sampling, 
based  on  the  sampling  pattern  and  on  the  time  intervals  between  the  hydrographic  stations. 

The  vertical  motion  was  also  diagnosed  from  the  analyzed  density  field  within  the 
framework  of  quasigeostrophic  (QG)  dynamics  by  solving  the  Q-vector  form  of  the  QG 
omega  equation.  The  results  thus  obtained  are  in  close  agreement  with  the  PE/DFI  results. 
The  same  patterns  of  rising  and  sinking  motion  were  again  found  and  explained  using  QG 
dynamics.  The  maximum  discrepancy  between  the  DFI  and  the  QG  vertical  motion  is  a 
magnitude  difference  of  about  30%,  with  the  QG  method  presenting  the  larger  values.  A 
clear  explanation  for  this  was  given,  based  on  the  terms  neglected  in  the  QG  omega 
equation  (the  two  most  important  in  this  study  are  thought  to  be  the  ageostrophic 
advections,  and  neglecting  relative  vorticity  in  the  divergence  term)  and  illustrated  using 
an  example  of  the  actual  ageostrophic  velocity  field  in  the  vicinity  of  one  of  the  meanders. 

From  all  of  the  results  summarized  above,  it  is  concluded  that  a  robust  description 
of  the  circulation  in  the  observed  cyclonic  eddy,  including  the  vertical  motion,  has  been 
achieved,  complementing  significantly  the  wealth  of  information  that  has  become  available 
with  recent  observational  programs  in  the  Californian  and  Oregonian  coasts.  The  detailed 
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circulation  in  eddies  such  as  the  one  studied  here,  i.e.,  with  the  same  temporal  and 
horizontal  scales,  and  at  the  same  stage  in  their  life  span,  can  be  inferred  from  the  results 
above.  The  methodology  used  in  this  work  can  also  be  applied  to  other  areas  of  the  world 
ocean  where  appropriate  quasisynoptic  data  sets  are  available. 

Future  work  should  be  directed  at  the  inclusion  of  bottom  topography  in  the  PE 
model  domain.  The  PE  results  in  the  region  of  the  domain  closer  to  shore  are  not  reliable 
because  the  flat  bottom  assumption  forced  some  of  the  profiles  to  be  extended  to  the 
uniform  depth  of  2000  m  using  a  similarity  method.  The  effect  of  a  varying  bottom  in  the 
vorticity  budget  will  surely  yield  a  more  complete  understanding  of  the  dynamics  of  the 
eddy,  and  eliminate  the  undesirable  effect  of  the  extended  profiles. 

Other  points  that  can  be  investigated  deal  with  the  optimization  of  the  DFI 
procedure,  namely  the  optimization  of  the  lowpass  filter  and  of  the  filtering  process,  since 
it  was  shown  that  filters  with  fewer  coefficients  can  better  meet  the  constraints  of  the  spin- 
up  phase. 

Another  aspect  that  is  worth  investigation  is  to  predict  the  future  motion  of  the  eddy . 
Once  the  fields  have  been  dynamically  balanced  by  the  DFI,  then  the  PE  model  can  be  used 
to  integrate  forward  and  simulate  the  evolution  of  the  eddy.  A  potential  problem  in  this 
particular  data  set  is  the  absence  of  information  on  the  temporal  evolution  of  the 
boundaries.  If  the  boundary  conditions  are  chosen  to  remain  constant  at  inflows  (as  they 
were  during  the  short  initialization  runs),  then  the  eddy  will  not  be  able  to  drift  westward 
and  the  general  result  will  perhaps  be  an  elongated  eddy  that  remains  attached  to  the 
boundaries  where  inflow  takes  place.  Thus  a  careful  treatment  of  the  boundaries  would  be 
necessary. 
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